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ABSTRACT 

This  thesis  deals  primarily  with  the  evaluation  of 
the  effects  of  the  changes  in  reservoir  water  levels  on 
the  stability  of  natural  slopes.  For  this  purpose  a  de¬ 
tailed  study  of  the  steady  state  unconfined  ground  water 
flow  through  a  porous  medium  is  undertaken. 

Finite  difference  techniques  have  been  used  to  deter¬ 
mine  the  phreatic  surface  and  potentials  in  a  slope  con¬ 
sisting  of  a  homogeneous,  isotropic  material  with  typical 
boundary  conditions. 

A  Fortran  IV  program  based  on  numerical  techniques 
has  been  written  to  solve  the  steady  state  free  surface 
flow  with  infiltration.  The  program  is  versatile  and  can 
be  used  for  any  geometric  configuration  and  boundary  condi¬ 
tions  with  or  without  infiltration.  This  program  has  been 
used  in  the  evaluation  of  the  amount  of  outflow  from  the 
banks  into  the  reservoir  at  various  reservoir  water  levels, 
in  addition  to  determining  the  pore  pressure  data  required 
for  checking  the  stability  of  slopes. 

Dimensionless  curves  have  been  developed  for  a  typical 
geometric  configuration  of  a  slope  which  can  be  used  for 
calculating  the  amount  of  outflow  from  the  banks  into  the 
reservoir  at  various  water  levels,  if  the  properties  of  the 
medium  are  known  and  if  the  initial  head  at  the  inside 
boundary  is  kept  constant. 

A  hypothetical  concept  of  constant  discharge  under 


steady  state  seepage  conditions  which  results  in  the  rais¬ 
ing  of  the  initial  head  at  the  inside  boundary  with  the 
raising  of  reservoir  water  level  is  postulated  and  graphs 
are  developed  to  determine  the  change  in  head  required  to 
keep  the  ground  water  discharge  constant  at  various  reser¬ 
voir  water  levels.  The  stability  of  slopes  is  then  analyzed 
and  results  are  compared  with  the  cases  of  constant  head. 

The  ICES  computer  program  LEASE-1  based  on  the 
simplified  Bishop  slip  circle  analysis  has  been  used  for 
checking  the  stability  of  slopes. 

The  influence  of  infiltration  on  the  phreatic  surface 
and  on  the  stability  of  slopes  is  examined.  Typical  re- 
sul ts  are  compared . 

Finally,  on  the  basis  of  the  constant  discnarge  concept, 
a  possible  explanation  of  Vajont  rock  slide  is  given. 
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CHAPTER  I 


INTRODUCTION 

1  . 1  Purpose  of  the  Research 

The  stability  of  natural  slopes  subjected  to  seepage 
forces  is  dependent,  among  other  factors,  on  the  pore 
pressures  within  the  soil  or  rock  mass.  The  determination 
of  these  pore  pressures  in  turn  depends  on  the  definition 
of  flow  domain. 

A  variety  of  slope  stability  problems  associated 
with  ground  water  flow  demand  a  complete  knowledge  of  the 
geohydro  1 ogi cal  charac teri s ti cs  of  the  area.  The  very 
occurrence  and  movement  of  the  ground  water  which  is  an 
important  phase  of  the  hydrological  cycle,  depends  mainly 
on  the  geological  formation  and  the  topography  of  the  area. 
Almost  all  ground  water  is  meteoric  water  derived  from 
precipitation  by  the  process  of  infiltration.  The  other 
main  sources  of  ground  water  recharge  are  percolation  from 
influent  streams,  lakes,  reservoirs  and  glaciers.  A  sche¬ 
matic  representation  of  the  hydrological  cycle,  occurrence 
and  movement  of  ground  water  is  illustrated  in  Fig.  1.1. 

In  a  natural  landscape  ground  water  moves  under  the 
action  of  gravity  from  higher  to  lower  elevations.  With¬ 
out  any  external  i n te rfe ren ce  ,  a  ground  water  basin  is 
filled  with  water.  The  excess  is  discharged  and  it  may 
emerge  as  surface  water.  Depending  upon  the  geohy drol ogi ca 1 
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FIG.  1.1  THE  HYDROLOGIC  CYCLE  (AFTER  TODD,  1959) 
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conditions  and  the  local  topography  the  discharge  of  ground 
water  tends  to  maintain  a  balance  between  the  inflow  and 
outflow.  It  is  therefore  reasonable  to  conclude  that  if 
the  seasonal  variations  of  the  natural  recharge  are  not 
very  significant,  an  aquifer  may  yield  a  constant  outflow. 
This  may  be  evident  as  a  head  water  source  of  a  small 
stream  or  merely  an  effluent  seepage  which  may  evaporate 
from  the  ground  surface. 

The  steady  state  unconfined  flow  pattern  in  a  region 
may  be  determined  by  seeking  the  solution  of  Laplace's 
equation.  This  is  done  by  using  one  of  the  available  tech¬ 
niques  with  the  prescribed  initial  and  boundary  conditions. 
The  general  practice  of  determining  the  free  surface  line 
is  to  fix  the  initial  head  well  inside  the  slope  and  the 
exit  conditions  at  the  outside  boundary.  This  initial  head 
which  is  usually  based  on  piezometric  observations  is 
chosen  at  a  location  well  inside  the  slope  and  is  assumed 
to  be  uninfluenced  by  any  change  in  the  geometric  configura¬ 
tion  of  the  slope  and/or  the  exit  conditions  of  the  flow. 
However,  in  order  to  maintain  a  balance  of  inflow  and  out¬ 
flow  which  is  an  inherent  cycle  of  recharge  and  discharge 
of  ground  water,  any  change  in  exit  conditions  of  flow  at 
the  outer  boundary  must  influence  the  initial  head  at  the 
inside  boundary  under  steady  state  seepage  conditions. 

The  impounding  of  a  reservoir  results  in  a  change  of 
the  exit  conditions  of  the  natural  ground  water  movement. 

If  the  initial  head  is  assumed  to  remain  constant  for  differ 
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ent  reservoir  water  levels,  the  phreatic  surface  would 
only  shift  with  respect  to  this  boundary  condition. 
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This  assumption  of  constant  head  leads  to  the  con¬ 
clusion  that  the  raising  of  the  reservoir  water  level  would 
result  in  the  reduction  of  the  outflow,  which  is  at  variance 
with  reality  for  most  of  the  field  cases  under  steady  state 
seepage  conditions. 

Transient  conditions  are  always  on  the  safe  side  as 
far  as  the  stability  of  slopes  is  concerned. 

It  therefore  follows  that  the  basic  assumption  of 
constant  head  does  not  hold  for  different  reservoir  water 
levels  under  steady  state  seepage  conditions,  whereas  a 
concept  of  constant  discharge  seems  to  be  more  realistic. 

In  order  to  assess  the  qualitative  and  quantitative 
implications  of  this  hypothesis  on  the  stability  of  natural 
slopes,  a  complete  study  of  the  steady  state  free  surface 
ground  water  flow  with  more  realistic  initial  and  boundary 
conditions,  defining  the  actual  flow  domain  in  the  field 
is  necessary. 
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1 . 2  Review  of  Previous  Work  Done 

The  very  occurrence  and  movement  of  ground  water 
depends  mainly  on  the  geohydrol ogi cal  conditions  of  the 
area  and  the  local  topography.  It  is  an  important  phase 
of  the  hydrological  cycle  and  without  any  external  inter¬ 
ference  maintains  a  balance  between  the  inflow  and  outflow 
(Todd  (  1959);  Linsley,  J.R.  e_t  aj_.  (  1958)). 

The  problem  of  unconfined  ground  water  flow  is  mostly 
treated  as  a  two  dimensional  problem.  Several  techniques 
exist  for  solving  the  free  surface  problems  of  steady  state 
unconfined  saturated  flow  through  porous  media. 

Earlier  work  is  based  on  graphical  methods.  These 
methods  of  flow-net  sketching  known  as  Dupuit's  solution 
(1863);  Schaffernak  and  Van  Iterson's  solution  (1917);  L. 
Casagrande's  solution  (1932);  and  Pavlosky's  solution  (1931) 
are  based  on  the  Dupuit  assumptions  and  give  approximate 
solutions. 

Some  analytical  solutions  as  reviewed  by  Harr  (1962), 
they  are  provided  by  Pol ubari nova-Kochi na  (1962);  Aravin 
and  Numerov  (  1  965);  and  Bear  e_t  al_.  (  1968).  All  these 
methods  become  too  complicated  to  be  applied  in  field  cases 
of  difficult  boundary  conditions  and  geometric  con f i gurati ons . 

Numerical  techniques  making  use  of  high  speed  digital 
computers  have  been  developed  to  overcome  these  difficul¬ 
ties.  Finite  difference  methods  using  the  relaxation  tech¬ 
niques  of  Allen  (1954);  Shaw  and  Southwell  (1941);  McNown 
et  al .  (  1  953)  and  Thom  and  Apel  t  (1961)  have  been  used  by 


various  investigators  in  solving  the  steady  state  free 
surface  flow  problems. 

Jeppson  (1968);  Finnemore  and  Perry  (1968)  have 
adopted  these  relaxation  techniques  in  analyzing  the  seep¬ 
age  through  an  earth  dam.  Jeppson  (1968)  and  Charmonman 
(1967)  have  also  applied  these  finite  difference  methods 
in  analyzing  the  seepage  through  canals  or  ditches. 

The  steady  state  free  surface  seepage  problems  have 
also  been  solved  by  Taylor  and  Brown  (1967);  Finn  (1967); 
Zienkiewicz  e_t  aj_.  (  1  966  );  Zienkiewicz,  0.  and  Cheung,  Y.K. 

(  1965  );  Neuman  and  Witherspoon  (  1  970)  and  Cheung  e t  a  1 . 
(1970)  by  the  finite  element  method. 

The  numerical  techniques  given  by  the  aboye  mentioned 
authors  using  either  finite  difference  or  finite  element 
methods  may  be  used  for  defining  the  free  surface  flow  with¬ 
out  any  infiltration  over  the  free  surface. 

Pol ubari nova-Kochi na  (1962)  has  given  an  exact  analy¬ 
tical  solution  for  analyzing  free  surface  seepage  with  con¬ 
stant  infiltration  over  the  entire  free  surface.  Verruijt 
(1970)  has  also  given  an  exact  solution  based  on  complex 
variable  techniques.  These  solutions  are  again  too  diffi¬ 
cult  to  apply  to  the  field  cases  of  complicated  boundaries 
and  geometric  configurations. 

The  stability  analyses  of  natural  slopes  subjected 
to  seepage  forces  requires  the  knowledge  of  the  geohydrology 
of  the  area,  the  properties  of  the  material  and  the  boundary 
conditions.  The  method  of  calculating  the  stability  depends 
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on  the  properties  of  the  material,  the  geometric  configura¬ 
tion  and  the  type  of  analysis,  i.e.,  short  term  or  long 
term  stability  analysis.  It  has  been  shown  by  Bishop 
(1952);  Henkel  and  Skempton  (1955)  and  Bishop  and  Bjerrum 
(1960)  that  the  effective  stress  analysis  is  more  advanta¬ 
geous  when  dealing  with  the  long  term  stability  of  slopes. 

A  sufficiently  accurate  method  based  on  effective 
stresses  is  given  by  Bishop  (1955),  who  treats  the  failure 
surface  as  an  arc  of  a  circle.  Although  the  Bishop  simpli¬ 
fied  method  does  not  satisfy  fully  statics,  it  has  been  shown 
that  application  of  the  Morgenstern  and  Price  (1965)  method 
to  circular  slip  surfaces  gives  app roxi matel y  the  same 
results.  If  the  material  is  homogeneous  and  isotropic, 
the  Bishop  simplified  method  may  be  used  without  any  appre¬ 
hension  . 

Almost  all  the  studies  dealing  with  the  stability  of  a 
slope  subjected  to  seepage  forces  adopt  the  general  practice 
of  fixing  the  initial  head  at  the  upstream  boundary  well 
inside  the  slope.  This  is  generally  based  on  piezometric 
observations  and  it  is  assumed  that  this  head  is  not  in¬ 
fluenced  by  any  change  in  the  exit  and/or  geometric  configura¬ 
tion  of  the  slope. 

As  an  example  it  is  of  interest  to  note  the  compre¬ 
hensive  work  on  the  rock  slide  in  the  Vajont  valley  which 
has  been  reported  by  various  investigators.  All  the  detailed 
studies  dealing  with  the  stability  of  the  Vajont  rock  slide 
after  this  unfortunate  disaster  are  more  or  less  based  on 


8 


the  assumption  of  constant  head  conditions. 

Muller  (1964,  1968)  in  describing  the  factors  influenc¬ 
ing  the  course  of  movement  had  considered  the  joint  water 
pressure  as  a  result  of  a  sloping  water  table  towards  the 
lake.  According  to  him,  in  1961,  the  mountain  ground  water 
level  sloped  at  an  inclination  of  about  4:1  towards  the 
side  of  the  valley,  correspond!' ng  to  a  hydrostatic  thrust 
of  the  joint  water  on  the  creeping  mass,  estimated  at  about 
2  to  4  million  tons.  This  estimate  is  however  based  on  the 
assumption  that  the  slope  of  the  ground  water  table  remains 
the  same  even  at  higher  lake  levels.  Later  on,  on  the 
basis  of  limited  piezometric  observations  it  was  concluded 
that  complete  levelling  of  the  inclination  of  the  ground 
water  occurred. 

Kenney  (1967)  in  calculating  the  maximum  unstabiliz¬ 
ing  effect  of  raising  the  reservoir  level  at  Vajont,  also 
assumed  that  the  water  level  in  the  rock  material  prior  to 
failure  was  horizontal  and  equal  to  the  reservoir  level. 

Jaeger  (1969)  has  discussed  the  stability  of  partly 
immersed  fissured  rock  masses  with  reference  to  the  Vajont 
rock  slide.  He  has  accepted  the  adverse  influence  of  the 
sloping  water  table  on  the  stability,  but  while  analyzing 
the  Vajont  rock  slide  has  assumed  that  the  water  level  in 
the  fissured  rock  is  about  the  same  as  in  the  steady  re¬ 
servoir. 

Since  the  overall  geohydrology  of  the  area  was  not 
known  and  the  piezometric  observations  were  very  limited 
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(e.g.,  only  three  piezometric  observation  locations  and 
the  farthest  one  was  located  at  an  elevation  less  than 
900  m.  and  approxi mate! y  800  m.  away  measured  horizontally 
from  the  deepest  point  of  the  river  (Muller  (1964)),  this 
assumption  of  levelling  of  water  table  was  made  without 
any  due  consideration  of  the  influence  of  the  filling  of 
the  reservoir  on  mountain  ground  water  table  further  from 
the  face  of  the  si  ope . 

During  the  first  stage  of  filling,  the  sliding  move¬ 
ments  started  at  the  reservoir  elevation  635  m.  and  became 
quite  high  (8-10  cms/day)  at  the  reservoir  elevation  650  m. 
As  pointed  out  by  Skempton  (1966),  the  reasons  for  this 
first  sliding  movements  have  to  be  explored.  These  move¬ 
ments  stopped  when  the  reservoir  level  was  lowered  to  600  m. 
and  remained  practically  negligible  during  the  period  of 
constant  reservoir  level  as  well  as  during  the  second  stage 
of  filling  up  to  the  reservoir  level  of  650  m.  The  movements 
were  magnified  again  when  the  previous  reservoir  level  was 
crossed  (see  Muller,  1964,  Fig.  19,  p.  174).  It  is  also 
evident  from  Muller,  1964,  Fig.  19,  p.  174  that  no  piezo¬ 
metric  data  of  the  area  well  inside  the  slope  are  available 
before  October,  1961.  Muller  (1968)  has  explained  this 
first  time  slide  by  considering  an  artesian  thrust  at  the 
base  of  the  sliding  mass,  which  is  based  on  the  assumption 
that  the  strata  lying  under  the  slip  surface  were  more  per¬ 
meable  than  the  strata  above. 

Jaeger  (1969)  has  explained  this  first  time  slide 


■ 


* 


10 


by  assuming  low  values  of  angle  of  shearing  resistance 
along  the  upper  strata  and  assuming  the  ground  water  table 
as  horizontal. 

Several  authors  have  calculated  the  angle  of  shearing 
resistance  at  failure.  Depending  upon  the  profile  chosen 
these  values  range  from  17.5°  to  24°  (Muller  (1968)). 
Nonveiller  (1967)  has  obtained  values  of  angle  of  shearing 
resistance  at  failure  from  27°  to  28°.  As  pointed  out  by 
Muller  (1968)  these  values  correspond  to  the  position  and 
shape  of  slip  surface  which  differ  very  much  from  nature. 

These  values  of  angle  of  shearing  resistance  obtained 
at  failure  are  much  smaller  than  what  could  be  expected 
for  limestone.  According  to  Skempton  (1966)  for  limestone 
a  friction  angle  of  about  30°  should  be  expected. 

To  summarize,  it  may  be  stated  that  all  these  excellent 
references  have  tried  to  explain  this  slide  of  very  com¬ 
plex  nature,  but  there  are  some  questions  which  still  re¬ 
main  unresol ved . 


1  .3  Scope  of  the  Study 


This  study  is  concerned  with  the  determination  of 
steady  state  unconfined  ground  water  flow  domains  within 
a  homogeneous,  isotropic  natural  slope  with  represen tati ve 
boundary  conditions.  The  boundary  conditions  are  defined 
with  special  reference  to  the  filling  of  reservoirs.  This 
results  in  a  change  in  flow  domain  and  in  turn,  a  change  in 
pore  pressures.  The  effect  of  this  change  in  pore  pressures 
at  various  reservoir  water  levels  on  the  stability  is  eva- 
1 uated . 

The  scope  of  the  study  also  includes  a  summary  of 
techniques  available  for  solving  free  surface  ground  water 
f 1 ow  p rob! ems . 

Finite  difference  techniques  for  the  numerical  solu¬ 
tion  of  Laplace's  equation  have  been  used  to  solve  the 
free  surface  flow  through  porous  media  with  constant  in¬ 
filtration. 

The  Fortran  IV  program  given  here  is  general  and  can 
also  be  used  for  determining  free  surface  flow  without  in¬ 
filtration.  Typical  results  computed  for  the  problem  of 
seepage  through  a  dam  without  infiltration  are  in  close 
agreement  with  Numerov's  analytical  solution  (Harr  (1962)) 
(see  Fig.  1.2). 

The  influence  of  raising  of  reservoir  water  level  on 
the  steady  state  seepage  and  consequently  on  the  stability 
of  slopes  is  studied  in  detail. 

The  constant  discharge  hypothesis  which  is  more  rea- 
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FIG.  1.2  COMPARISON  OF  ANALYTICAL  AND  NUMERICAL  SOLUTIONS 
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listic  than  constant  head  for  field  cases,  is  postulated 
and  studies  have  been  made  to  assess  the  influence  of  this 
on  the  stability  of  slopes. 

A  study  of  outflow  from  the  banks  into  the  reservoir 
for  various  reservoir  water  levels  is  also  made  and  dimen¬ 
sionless  curves  are  presented,  which  can  be  used  for  cal¬ 
culating  the  seepage  through  dams  with  different  tail  water 
levels  and/or  outflow  from  the  banks  into  the  reservoir, 
knowing  the  properties  of  soil  and  keeping  the  initial 
boundary  head  as  constant. 

The  influence  of  infiltration  on  the  free  surface 
flow  and  then  on  the  stability  of  slopes  is  examined  and 
results  are  compared  with  those  of  steady  seepage  conditions 
without  infiltration. 


' 


CHAPTER  II 


HYDRAULIC  CONCEPT  OF  FREE  SURFACE 


GROUND  WATER  FLOW 


2 . 1  General 

The  problem  of  unconfined  ground  water  flow  will  be 
treated  in  two  dimensions. 

Assuming  the  seepage  flow  to  be  laminar,  which  is 
usually  the  case  in  practice,  the  applicability  of  Darcy's 
law  is  valid  and  the  following  fundamental  equations  for 
two  dimensional  flow  in  the  xy  plane  are  obtained. 


(2.1) 


(2.2) 


where  u  and  v  are  the  components  of  discharge  velocity  in  the 
x  and  y  direction  respecti vely . 

Taking  the  y  axis  as  vertical,  oriented  upward,  and 
the  x  axis  horizontal,  the  head  h  and  velocity  potential 
<J>  may  be  defined  as 


h 


(2.3) 
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<i>  =  -k(z  +  £-)  +  c  (2.4) 

Yw 

where  z  is  the  elevation  head  of  the  point  above  an  arbi¬ 
trary  reference  datum,  P  is  the  pressure  in  the  fluid  at 

that  point,  k  is  the  coefficient  of  permeability  and  y,f  is 

w 

the  specific  weight  of  water. 

The  continuity  equation  may  be  expressed  as 


3  u 

3x 


+ 


(2.5) 


Replacing  u  and  v  in  the  equation  of  continuity,  the  well 
known  Laplace  equation  is  obtained. 


(2.6) 


Assuming  the  potential  flow  to  obey  Laplace's  equa¬ 
tion,  the  lines  of  equal  velocity  potentials  are  orthogonal  to 
the  flow  lines.  In  ground  water  flow  another  harmonic  function 

known  as  the  stream  function  ip  ( x ,  y )  is  introduced,  which 
is  defined  as 


v 


3x 


(2.7) 
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The  stream  function  ^(x,y)  satisfies  Laplace's  equation 
if  it  is  defined  by  the  Cauchy- Ri emann  equations.  The 
Cau chy- Ri emann  equations  are  obtained  by  combining  equa¬ 
tions  (2.1),  (2.2)  and  (2.7). 


8 d>  _  dip 
8x  8y 

(2.8) 

d<p  _  dip 
dy  ~  dx 


Since  by  definition  it  follows  that  the  velocity 
potential  is  a  single  valued  function  at  every  point  of 
the  x-y  plane,  it  is  possible  to  draw  lines  of  constant  <j> 
in  the  x-y  plane.  These  lines  are  called  potential  lines 
and  are  generally  drawn  at  constant  intervals  (Fig.  2.1). 


(J)-j  — 1  (j)  ^  $3  ...  Ac|) 

t 

Referring  to  Fig.  2.1  the  directions  s  and  n  are 
defined  by 


s  =  x  cos  a  +  y  sin  a 


n  =  -  x  sin  a  +  y  cos  a 


(2.9) 


Simplifying  and  rearranging,  we  get 


' 
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FIG.  2.1 


representation  of  potential  and  flow  lines 
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x  =  s  cos  a  -  n  sin  a 


(2.10) 


y  =  s  sin  a+  n  cos  a 


The  angle  a  is  chosen  in  such  a  w ay  that  the  n-axis  is 
tangent  to  the  potential  line  and  the  s-axis  is  perpendi¬ 
cular  to  it.  It  follows  from  the  tangency  of  the  n-axis 
with  the  potential  line  that 


(2.11) 


Hence  there  is  no  specific  discharge  component  in  the 
n-direction.  This  obviously  means  that  the  direction  of 
flow  is  perpendicular  to  the  potential  lines  and  flow  occurs 
only  in  the  s-direction. 

If  the  magnitude  of  the  specific  discharge  vector  is 
denoted  by  q  ,  we  have 


q 


(2.12) 


The  velocity  components  expressed  in  q  are 


q  =  q  cos  a 

X 


(2.13) 


qy  =  q  sin  a 
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Darcy's  law  may  also  be  written  as 


(2.14) 


and  the  continuity  equation  in  terms  of  q  becomes 


!!x 

dx 


+ 


(2.15) 


It  follows  from  the  continuity  equation  that  the  flow 
can  also  be  derived  from  the  single  valued  stream  function 
i l>(x,y)  or  simply  \p  by 


q 


x 


(2.16) 


Since  equation  (2.15)  is  satisfied  identically  it  follows 
from  equation  (2.14)  that 


ay  3x 


(2.17) 
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Substitution  of  equati  on  (2 . 1 6 )  in  (2.17)  gives 

V2iJ>  =  (2.18) 

3x  3y 

which  shows  that  ip  is  also  a  harmonic  function. 

It  is  interesting  to  note  that  4>  does  not  vary  in  the 
s-direction,  hence  the  direction  of  flow  is  tangent  to  the 
lines  of  constant  ip.  This  can  easily  be  demonstrated  as 
f ol 1 ows  : 


ii  =  M.  dx  +  3i  dy_ 
9s  9x  ds  9y  ds 


=  q  cos  a  -  qx  sin  a 

V 


(2.19) 


=  q  sin  a  cos  a  -  q  cos  a  sin  a 


=  0  . 


Si  mi  1 arl y 


3 \b  _  djp  dx 
9n  9x  dn 


11  ly 

9y  dn 


=  -  q  sin  a  -  q  cos  a 

y 

(2.20) 

2  2 
=  -  q  sin  a  -  q  cos  a 

=  -  q. 
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A  comparison  of  equation  (2.12)  and  (2.20)  shows 

that 


djp  _  dip 
3s  3n 


(2.21) 


This  means  that  in  the  x~y  plane  if  lines  of  constant  cp  and 
ip  are  drawn  at  intervals  A <p  and  At/;,  these  lines  will  form 
squares  and  if  we  choose 


then 


As  =  An 


A(J>  =  Aip  . 

The  total  discharge  AQ  between  two  stream  lines  at  a  mutual 
distance  of  An  may  be  given  by 

AQ  =  qBAn  =  -  B  A  if;  (2.22) 

where  B  is  the  thickness  of  the  plane  where  flow  occurs. 

This  relationship  shows  that  the  discharge  per  unit 
thickness  between  two  stream  lines  equals  the  difference 
in  the  values  of  \p  of  these  two  stream  lines  and  may  be 
used  in  making  outflow  studies. 
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2 . 2  Mathematical  Formulation  of  the  Problem 

The  analysis  of  steady  flow  through  porous  media, 
involving  the  determination  of  the  flow  pattern  and  the 
potential  distribution  re qu ires  the  solution  of  Laplace's 
equation  subject  to  representa ti ve  boundary  conditions. 

This  normally  necessitates  a  detailed  knowledge  of  ground 
water  flow  characteristics,  the  boundary  conditions  and  the 
change  in  the  flow  pattern  due  to  any  changes  in  the  initial 
or  boundary  conditions.  It  is  therefore  essential  to  define 
explicitly  the  initial  and  boundary  conditions  representa¬ 
tive  for  the  field  case  under  study. 

A  hypothetical  case  of  a  natural  slope  consisting  of 
homogeneous  and  isotropic  material,  subjected  to  various 
water  pressures  on  the  down  stream  side,  corresponding  to 
the  filling  of  the  reservoir  is  studied  in  detail  (Fig.  2.2a). 
The  following  boundary  assumptions  are  made: 

(1)  A  known  potential  at  the  upstream  boundary.  This  up¬ 
stream  boundary  is  generally  selected  at  a  location 
where  the  total  head  is  assumed  to  be  constant.  This 
means  that  along  the  boundary  AE  the  potential  lines  are 
vertical  and  the  pressure  distribution  is  hydrostatic. 

(2)  A  known  potential  at  the  down  stream  boundary  CD,  which 
is  also  a  line  of  constant  potential. 

(3)  An  impervious  boundary  ED  where  no  flow  takes  place 
across  it  and  thus  ED  is  a  locus  of  a  stream  line. 

(4)  Free  water  table  AC,  which  is  the  upper  limit  of  the 
seepage  in  the  flow  domain.  Along  this  line  the  pres- 
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FIG.  2.2a  BASIC  LAYOUT  OF  SLOPE 
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sure  at  every  point  is  equal  to  atmospheric  pressure 
and  that  this  free  surface  is  a  stream  line. 

(5)  When  dealing  with  a  case  of  constant  infiltration  over 
the  entire  free  surface  some  new  boundary  conditions 
have  to  be  introduced.  A  detailed  account  of  the 
typical  boundary  conditions  and  the  derivation  of  the 
resulting  governing  differential  equation  for  the  free 
surface  with  infiltration  is  given  as  follows: 

Assume  a  free  surface  line  sloping  at  an  angle  a  with 
the  horizontal  and  having  a  constant  infiltration  over  the 
entire  surface  (Fig.  2.2b(i) and  Fig.  2.2b(ii ))  .  If  q  is 
the  volume  of  water  reaching  the  free  surface  per  unit  area 
in  the  horizontal  plane,  per  unit  time,  we  have 


M  = 

3n 


-  q  COS  a 


(2.23) 


where  a  is  the  angle  measured  clockwise  from  the  horizontal 
to  the  tangent  of  the  free  surface.  When  a  =  0,  i.e.,  the 
free  surface  line  is  horizontal 


3 

3n 


-  q 


(2.24) 


when  a  =  90°,  i.e.,  the  free  surface  line  is  vertical 


(2.25) 


. 
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FIG.  2. 2b( i )  FREE  SURFACE  WITHOUT  INFILTRATION 


FIG.  2.2b( i i )  FREE  SURFACE  WITH  INFILTRATION 
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we  also  know  that 


3(j>  _  3tfr  dx  +  3<j>  dy 

3n  9x  dn  3y  dn 


3d)  •  .  3  (b 

J  sin  a  +  J-  cos  « 


It  follows  from  equation  (2.23)  that 


3d)  •  ,3  d) 

—  s l n  a  +  cos  a  = 


-  q  cos  a 


(2.26) 


we  also  know  that  at  every  point  along  the  free  surface 
line  the  pressure  should  be  atmospheric,  so  that 

(p  +  ky  =  constant  . 

Differentiating  along  s,  we  have 


|±  +  k  &  =  o 

3s  ds 


(2.27) 


3<j)  3(f)  dx  +  3(f)  dy 

3s  3x  ds  3y  ds 


If  cos  a  -  f-f  sin  a 


Again,  we  know 
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Substituting  in  equation  (2.27),  we  get 


3d)  3d)  , 

—■cos  a  -  yj-  s  i  n  a  =  k  sin  a 


(2.28) 


Rearranging  equations  (2.26)  and  (2.28) 


3d>  . 

3?  Sin  a  = 


-  (|y  +  q)  COS  a 


(2.29) 


and 


cos  a  =  (|y  +  k)  sin  a 


(2.30) 


It  follows  from  equation  (2.29)  and  (2.30) 


tan  a  = 


_  (M  +  q)  /II 

^  3y  Mi,/3x 


(2.31  ) 


and  also 


tan  a  = 


3xM3y  ' 


Eliminating  a  from  equation  (2.29),  we  get 


-  (ff)2  =  k  |i  +  q  +  (|£)2  +  qk 


3y  11  ay  'ay 


or 
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(ff)2  +  (|y)2  +  (k  +  q)  +  qk  =  0  (2.32) 

In  the  case  of  no  infiltration  equation  (2.32)  reduces  to 


(2.33) 


Equation  (2.32)  and  (2.33),  which  apply  along  the  free  sur¬ 
face,  are  the  governing  differential  equations  for  the  free 
surface  with  or  without  infiltration  respectively  and  are 
identical  to  the  equations  given  by  Pol ubari nova-Kochi na 
(  1  962)  . 

In  terms  of  total  head  equation  (2.26)  may  be  written 
as 


9H 

3x 


s  i  n 


a  + 


3H 

3y 


cos 


a 


Q  cos  a 


(2.34) 


where  Q  is  defined  as  the  ratio  of  the  volume  of  water  in¬ 
filtrating  per  unit  horizontal  area  per  unit  time  divided 
by  the  coefficient  of  permeability  of  the  soil,  i.e., 

Q  =  -  q/k. 
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2 . 3  Numerical  Approach 

Most  of  the  free  surface  ground  water  problems  with 
complicated  geometry  and/or  boundary  conditions  become  too 
cumbersome  to  be  handled  analytically.  An  alternative  to 
analytic  methods  is  the  numerical  approach,  which  reduces 
the  problem  to  the  solution  of  a  set  of  algebraic  equations. 

It  has  already  been  demonstrated  by  Southwell  (1946), 
Shaw  (1953),  McNown  e_t  aj_.  (1953)  and  Allen  (1954)  that 
free  surface  problems  can  be  handled  by  finite  difference 
methods.  Various  schemes  are  available  to  transform  the 
differential  equations  into  a  finite  difference  form.  The 
Laplace  equation  is  generally  approximated  by  the  following 
finite  difference  equation 

H  ( I  ,  J )  =  (H(I,J  +  1)  +  H(I-1,J)  +  H  ( I  ,  J  - 1  ) 

+  H ( I  +  l  ,J))  *  0.25  (2.35) 


where  H  (I,J)  represents  the  total  head  at  the  grid  point 
( I , J )  (Fig.  2.3).  This  finite  difference  equation  has  been 
obtained  by  using  a  five  point  grid  system,  a  square  network, 
and  is  therefore  based  on  a  second  degree  polynomial  approxi¬ 
mation.  It  is,  however,  possible  to  develop  finite  differ¬ 
ence  equations  with  smaller  truncation  errors  Thom  and  Apelt 
(1961),  but  in  view  of  the  complicated  algebra  involved  in 
combining  the  boundary  conditions,  the  most  promising  approach 
seems  to  be  a  five  point  grid  system. 
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ASSUMED  FREE  SURFACE  LINE 


FIG  2  3.  REPRESENTATION  OF  COMPUTATIONAL  MOLECULE  AND  FREE  SURFACE 
LINE  WITH  IRREGULAR  BOUNDARY 


The  finite  difference  equations  can  be  solved  by 
several  methods  including  direct  matrix  inversion,  the 
Gauss-Seidel  type  iterative  procedure,  etc.  By  these 
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methods  the  lengthy  "relaxation"  solution  is  obtained  re¬ 
peatedly  . 

Since  the  solution  of  these  equations  requires  an 
iterative  process,  the  most  efficient  method  is  the  Gauss- 
Seidel  iterative  method  with  an  o ve r-re 1 axa ti on  factor. 

This  method  ensures  rapid  convergence  and  thereby  less  num¬ 
ber  of  iterations  are  required.  This  o ve r- rel axati on  method 
is  known  as  the  "successive  over-rel axati on "  (SOR)  method 
and  can  be  expressed  mathematically  as 

Hk(I,J)  =  Hk~ 1  (I  ,0 )  +  |  *  [H  ( 1 ,0+1 )  +  H(I-l.J) 

+  H(I.J-l)  +  H(  1  +  1  ,H)  ]  -  M*Hk"1(I,J)  (2.36) 

in  which  k  is  the  kth  iteration  and  to  is  the  over- re  Taxation 
factor.  It  has  been  found  theoretically  that  a  value  of  w 
ranging  between  0  and  2  ensures  convergence  of  the  succes¬ 
sive  over-relaxation  method.  There  is,  however,  an  optimum 
value  of  a)  which  will  give  the  most  rapid  convergence 
(Forsythe  and  Wasow  (I960)). 

The  biggest  advantage  of  the  successive  over-relaxation 
method  is  that  it  does  not  require  any  storage  of  the  equa¬ 
tions,  but  the  value  of  the  function  is  adjusted  at  every 
grid  point  by  the  finite  difference  equation.  The  new  value 
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is  used  for  all  computations,  involving  the  values  obtained 
previously. 

In  order  to  determine  potentials  at  the  free  surface, 
second  degree  polynomials  were  tried.  These  polynomials 
approximate  the  function  between  three  consecutive  points 
in  each  local  interval.  These  points  need  not  lie  on  the 
grid  points  shown  on  Fig.  2.3;  in  that  case  linear  inter¬ 
polation  has  been  performed  to  obtain  the  values  of  poten¬ 
tials  at  these  poi nts . 

The  polynomial  tried  is  of  the  form 


H( x,y)  =  a  +  bx(X  -  XQ)  +  Cx(X  -  XQ)2 


+  VY  -  V  +  Cy(Y  "  Y0)2 


(2.37) 


It  follows  from  equation  (2.37)  that  (with  reference  to 
Fig.  2.3) 


(2.40) 


(2.38) 


(2.39) 


(2.41) 
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where  X-j  ,  and  Y  ,  correspond  to  the  abscissa  and 

ordinates  of  the  points  1,  2  and  a,  b  and  Xq  ,  Yq  are  the 
coordinates  of  the  point  0  on  the  free  surface  (Fig.  2.3). 

Vie  now  have  four  equations,  i.e.,  (2.38),  (2.39), 
(2.40)  and  (2.41)  and  five  unknowns.  Another  equation  is 
therefore  required  which  is  obtained  from  the  boundary 
condition  of  the  free  surface. 

Differentiating  equation  (2.37)  with  respect  to  x  and 
y  respecti vely  ,  we  get, 


=  bx  +  2Cx(X  -  x0)  (2.42) 

and 

fy  =  by  +  2Cy(Y  -  Y0)  (2.43) 

It  follows  from  equations  (2.42)  and  (2.43)  that  at  point 
0  on  the  free  surface 


9H 


Bx'O 


)n  =  b 


(2.44) 


3H 
9y '0 


)n  =  b 


(2.45) 


Therefore  equation  (2.34)  may  be  written  as 


b  sin  a  +  b,  cos  a  =  -  Q  cos  a 
x  y 


(2.46) 


The  angle  a  at  each  point  may  be  defined  as 


We  now  have  five  equations,  e.g.,  (2.38),  (2.39),  (2.40), 
(2.41)  and  (2.46)  which  are  set  up  in  matrix  form  and 
solved  for  <j>g  by  the  Gaussian  elimination  method. 

It  may  be  seen  that  equation  (2.34)  can  only  be  used 
for  interior  points  when  all  the  points  lie  on  the  grid 
point.  However,  the  free  surface  points  may  not  necessarily 
lie  on  the  grid  points  (Fig.  2.3).  In  that  case  different 
finite  difference  equations  based  on  a  five  point  grid 
system  with  irregular  boundaries  are  used.  There  may  be 
three  different  cases  corresponding  to  this  irregular  bound¬ 
ary  (Fig.  2.3)  which  are  described  below: 

(1)  When  the  free  surface  point  does  not  lie  on  the  grid 
point  but  cuts  the  vertical  grid  line  at  a  distance 
nA  above  the  center  grid  point  H(I,J)  and  all  the  other 
three  points  lie  on  the  grid  points,  the  following 
finite  difference  equation  is  used: 


Hc(  I  » J  )  [  n  (  n  +  ry 


2HS(J,1) 


2H(  1  +  1  ,J) 

1+n 


(2.47) 


+  H  ( I  ,0  +  1  )  +  H(  I  ,  J-l  )]  *  2(-”-+Tj 
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in  which  HS(J,1)  represents  the  potential  at  the  res¬ 
pective  points  on  the  free  surface  (Allen  (1954)). 

(2)  When  the  free  surface  point  cuts  the  horizontal  grid 
line  at  a  distance  AA  from  the  center  grid  point 
H(I,J)  and  the  other  three  grid  points  are  at  a  dis¬ 
tance  A  each  from  the  point  H(I,J),  the  finite  differ¬ 
ence  equation  used  is  of  the  following  form 


(2.48) 


+  H ( I  -1  ,0)  +  H(  1  +  1  ,J)]  *  2TI+T7 


(3)  Lastly  when  the  free  surface  line  cuts  the  vertical 

as  well  as  horizontal  grid  lines  at  a  distance  nA  and 
AA  respectively  from  the  center  grid  point,  the  follow¬ 
ing  finite  difference  equation  is  used 


2M  J+l  ,1 )  2 H o  ( J  ,  1  ) 

M1’^  =  r  -xiv+Tf—  +  w 


(2.49) 


2H ( I  , J -1  )  ,  2H( 1  +  1  ,J)  *  nA 

'  (1+A)  1+n  'J  2(n+A ) 


The  problem  of  the  impervious  boundary  is  treated  on 
the  basis  of  the  fact  that  there  is  no  flow  across  the 
boundary.  The  finite  difference  equation  obtained  for 
this  case  is  of  the  following  form: 
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hlc^IMAX’J)  "  rH(rMAX’J  +  1>  +  2H(!max-1’J) 


+  “(^ax.J-I)]  *  °-25 


(2.50) 
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2 . 4  Method  of  Solution 

The  correct  position  of  the  free  surface  is  obtained 
systematically  by  seeking  the  solution  of  the  finite  differ¬ 
ence  equations  correspondi ng  to  the  various  boundary  condi¬ 
tions. 

A  Fortran  IV  program  has  been  written  to  solve  this 
steady  state  free  surface  ground  water  flow  problem.  This 
program  determines  the  correct  position  of  the  free  sur¬ 
face  in  addition  to  the  potentials  at  the  grid  points. 

A  systematic  scheme  has  been  postulated  which  is  out¬ 
lined  stepwise  as  follows. 

$  tep  1 

The  location  of  the  free  surface  is  estimated  and  its 
coordinates  are  defined  on  the  grid  lines.  Initial  values 
of  potentials  at  the  free  surface,  at  the  interior  grid 
points  and  also  at  the  prescribed  upstream  and  downstream 
boundaries  are  given. 

$  tep  2 

Laplace's  equation,  approximated  by  finite  difference 
equations  is  solved  by  successive  over-relaxation  subject 
to  the  boundary  conditions.  The  residue  at  every  center 
grid  point  of  the  five  point  grid  system  is  determined  and 
checked  with  the  maximum  allowable  value.  The  iterative  pro¬ 
cess  is  con ti n ued  until  this  condition  is  satisfied. 
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S tep  3 

The  values  of  potentials  for  the  two  consecutive 
points  near  the  free  surface  which  are  not  lying  on  the 
grid  points  in  either  the  X  or  Y  directions  are  interpo¬ 
lated  and  are  used  in  solving  equations  (2.38),  (2.39), 
(2.40)  ,  (2.41 )  and  (2.46) . 

The  five  simultaneous  linear  equations  are  solved 
by  Gaussian  elimination  method  and  the  potentials  at  points 
on  the  free  surface  are  obtained. 

Step  4 

Using  these  values  of  potential  at  the  free  surface, 
the  potentials  at  the  interior  points  are  then  corrected 
by  using  SOR. 

The  process  of  relaxation  is  continued  till  the  speci¬ 
fied  value  of  the  residue  is  achieved.  Again  the  potentials 
at  free  surface  are  determined  in  a  similar  way  as  des¬ 
cribed  above. 

These  potentials  at  the  surface  are  compared  every 
time  with  those  obtained  from  the  previous  iteration  and 
the  whole  operation  is  repeated  until  the  difference  in 
potentials  at  all  the  free  surface  points  is  equal  to  or 
less  than  the  specified  value. 

S  te  p  5 

The  last  step  is  to  check  the  boundary  condition  at 
the  free  surface.  The  pressure  at  the  free  surface  should 
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be  atmospheric,  or  taken  to  be  equal  to  zero. 
We  know 


h  =  (—  +  y) 

yW 

or 


P  =  YW(H  -  Y) 

where  H  =  total  head  or  piezometric  head 
Y  =  elevation  head. 

For  the  condition  P  =  0,  H  must  be  equal  to  Y.  This  condi¬ 
tion  is  checked  at  every  point  on  the  free  surface  and  if 
it  is  not  satisfied,  the  value  of  Y  is  changed  at  the  res¬ 
pective  points.  This  means  a  new  trial  free  surface  is 
drawn  and  the  potentials  at  the  free  surface  are  determined. 
If  the  values  of  these  potentials  thus  obtained  at  the  res¬ 
pective  points  on  the  free  surface  are  equal  to  the  values 
of  Y  initially  given,  the  problem  is  solved  and  this  line 
represents  the  correct  free  surface. 

This  method  which  appears  to  be  lengthy  is  quite 
rapid.  The  trial  free  surface  lines  converge  very  rapidly 
and  the  final  solution  is  obtained  in  three  or  four  trials. 
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CHAPTER  III 


OUTFLOW  STUDY  AND  STABILITY  ANALYSES 

3 . 1  General 

The  amount  of  outflow  from  a  natural  slope  is  depen¬ 
dent  on,  among  other  things,  the  permeability,  the  slope 
of  the  ground  water  table  and  the  boundary  conditions.  The 
most  common  practice  is  to  fix  the  upstream  head  and  bound¬ 
ary  on  the  basis  of  some  piezometric  observations  prior  to 
the  filling  of  the  reservoir.  This  upstream  boundary  condi¬ 
tion  is  generally  kept  unchanged  in  determining  the  flow 
pattern  corresponding  to  various  reservoir  water  levels. 

If  the  slope  of  the  initial  equilibrium  water  table  is  not 
very  steep  and  the  upstream  boundary  is  chosen  at  a  location 
well  inside  the  slope,  it  might  be  a  fair  approximation  to 
assume  that  the  head  H  at  the  upstream  is  not  influenced 
by  the  changes  in  the  reservoir  water  levels.  On  the  other 
hand  if  the  geometric  con f i gura ti on  of  the  slope  and  the 
geohyd ro 1 o gi ca 1  conditions  of  the  area  are  such  that  the 
initial  slope  of  the  natural  ground  water  table  is  quite 
steep,  the  assumption  of  constant  head  for  determining  the 
pore  pressures  at  various  reservoir  water  levels  does  not 
seem  to  be  realistic. 

Various  methods  are  available  for  analyzing  the  stabi¬ 
lity  of  slopes.  Studies  have  shown  that  effective  stress 
analyses  are  more  useful  in  analyzing  the  long  term  stability 


. 
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of  slopes  and  embankments  (Bishop,  1952;  Henkel  and  Skempton, 
1955;  Bishop  and  Bjerrum,  1960).  It  therefore  follows  that 
in  order  to  use  this  analysis,  it  is  necessary  to  know  the 
pore  pressure  distribution  within  the  soil  or  rock  mass. 

A  sufficiently  accurate  method  of  stability  analysis 
based  on  effective  stresses  is  given  by  Bishop  (1955)  and 
is  known  as  Bishop's  simplified  method.  This  method  treats 
the  failure  surfaces  as  arcs  of  circles.  It  has  been 
adopted  in  this  study. 
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3 . 2  Outflow  Study 

The  main  aim  of  this  part  of  the  study  is  to  eva¬ 
luate  the  reduction  in  discharge  from  the  bank  with  the 
raising  of  reservoir  water  level  if  the  head  at  the  upstream 
boundary  is  kept  constant.  The  case  of  impounding  behind 
a  dam,  where  the  reservoir  water  level  is  kept  constant  is 
analogous  to  it.  Since  the  influence  of  the  upstream  slope 
of  a  dam  is  small  on  the  seepage  discharge  (Harr  (1962)), 
for  practical  purposes,  this  study  may  also  be  found  use¬ 
ful  in  calculating  seepage  through  a  dam  at  yarious  tail 
water  levels. 

For  this  purpose  a  hypothetical  case  of  a  30°  slope 
with  an  impervious  base,  consisting  of  a  homogeneous  and 
isotropic  material  has  been  adopted.  The  upstream  head 
is  fixed  at  an  elevation  of  1000  and  the  reservoir  water 
level  is  raised  from  0  to  800.  The  upstream  boundary  is 
chosen  at  a  distance  four  times  the  head  inside  the  bank 
and  is  assumed  to  be  uninfluenced  by  the  changes  in  the 
reservoir  water  levels.  The  computer  program  has  been  used 
to  determine  the  free  surface  and  potentials  for  each  case 
of  reservoir  water  level,  ranging  from  0  to  800  elevation 
in  increments  of  200  (Fig.  3.1  through  3.5).  Having  ob¬ 
tained  the  free  surface  line  and  the  potentials,  the  equi- 
potential  and  flow  lines  are  drawn  and  the  amount  of  outflow 
in  a  dimensionless  form  is  calculated  for  each  case.  A 
sample  calculation  of  q  for  the  case  d/H  =  0  is  given 

II I  cl  A 

in  Appendix  B. 
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A  plot  between  q/q  and  d/H  is  then  developed  and 

rn  a  x 

is  shown  on  Fig.  3.6.  In  this  plot  q  represents  the 

in  a  x 

discharge  when  there  is  no  water  in  the  reservoir,  q  is 
the  discharge  corresponding  to  the  respective  water  level 
d  in  the  reservoir  and  H  is  the  head  at  the  upstream  bound¬ 
ary.  As  evident  from  this  curve  the  relationship  between 
the  discharge  and  the  reservoir  water  levels  is  highly  non¬ 
linear.  The  reduction  in  discharge  is  not  very  significant 
at  low  reservoir  water  levels,  while  it  becomes  quite 
appreciable  at  higher  levels,  e.g.,  at  d/h  =  0.8,  the  re¬ 
duction  in  discharge  is  almost  50%. 
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q  =  Discharge  corresponding 
max  to  d=0 

q  =  Discharge  at  particular 
reservoir  level  d 
H  =  Head  at  the  U/S 
boundary 

d  =  Reservoir  level 


FIG.  3.6  PLOT  OF  DISCHARGE  AS  A  FUNCTION  OF  RESERVOIR  LEVEL 
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3 . 3  The  Constant  Discharge  Concept 

It  is  now  well  known  that  if  the  main  sources  of 
ground  water  recharge  are  precipitation,  glaciers,  influent 
streams,  lakes  and  reservoirs.,  there  exists  a  balance  bet¬ 
ween  inflow  and  outflow  in  the  ground  water  basin  under 
equilibrium  conditions.  This  results  in  a  certain  amount 
of  outflow  from  the  banks,  especially  in  hilly  terrain. 

This  outflow  may  vary  seasonally  but  if  all  the  sources  of 
recharge  have  a  cumulative  effect  on  the  occurrence  and 
movement  of  ground  water,  the  seasonal  variations  may  not 
be  very  significant.  However,  there  could  be  some  cases 
when  the  seasonal  variation  of  ground  water  table  is  quite 
appreciable.  In  that  case  the  safest  approach  for  stability 
analysis  is  to  consider  the  highest  and  steepest  ground 
water  table.  This  would  give  the  maximum  outflow  and  ob¬ 
viously  be  the  worst  case. 

In  most  of  the  field  cases  the  local  topography  plays 
an  important  role  in  controlling  the  regional  pattern  of 
the  ground  water  flow.  As  the  filling  of  the  reservoir 
proceeds  the  flow  pattern  is  altered.  Although  the  tran¬ 
sient  conditions  will  prevail  for  sometime,  the  final  steady 
state  conditions  are  more  realistic  for  all  but  the  more 
impermeable  materials.  If  the  head  at  the  upstream  boundary 
is  assumed  to  be  constant,  it  is  obvious  that  the  raising 
of  the  reservoir  water  level  would  result  in  the  reduction 
of  discharge  from  the  banks  into  the  reservoir.  This  con¬ 
clusion  does  not  seem  to  be  reasonable  as  according  to  the 
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concept  of  a  recharge  and  discharge  cycle  and  the  conserva¬ 
tion  of  mass,  the  amount  of  water  entering  at  the  upstream 
boundary  under  normal  conditions  has  to  go  out.  Further, 
if  there  is  a  continuity  of  flow,  the  raising  of  reservoir 
water  level  which  is  changing  the  exit  conditions  of  the 
ground  water  flow,  induces  a  raising  of  the  ground  water 
level  on  the  upstream  side.  This  will  ultimately  become 
asymtotic  to  the  original  water  table  at  infinite  distance. 
This  means  that  impounding  results  in  the  reduction  of  dis¬ 
charge  temporarily,  corresponding  to  the  transient  flow 
conditions,  but  under  steady  state  seepage  when  the  ground 
water  table  is  stabilized  again,  the  amount  of  water  flow¬ 
ing  out  should  remain  the  same.  Since  the  location  of  the 
upstream  boundary  is  not  changed,  the  permeability  of  the 
material  is  also  assumed  to  be  constant,  the  only  factor 
controlling  the  amount  of  outflow  at  various  reservoir 
water  levels  is  the  head  at  the  upstream  boundary.  Hence 
the  assumption  of  constant  head  is  not  valid  for  most  of 
the  field  cases;  on  the  contrary,  the  concept  of  constant 
discharge  seems  to  be  more  reasonable. 

To  evaluate  the  effects  of  this  hypothesis,  the  same 
problem  of  a  hypothetical  slope  with  typical  boundary 
conditions  as  explained  in  Section  3.2  is  studied.  The  free 
surface  line  and  the  potentials  for  the  case  when  the  re¬ 
servoir  is  empty  is  determined  (Fig.  3.1)  and  by  sketching 
flow nets  the  amount  of  outflow  from  the  slope  into  the 
reservoir  is  calculated  (see  Appendix  B). 
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The  increase  in  head  required  at  the  upstream  bound¬ 
ary  to  keep  the  discharge  constant  at  various  reservoir 
water  levels  is  then  found.  No  simple  relationship  could 
be  established  due  to  the  non-linear  relationship  between 
the  re  duct ion  in  discharge  and  the  increase  in  reservoir 
water  level.  As  an  alternative  the  method  of  trial  and 
error  is  followed  and  the  head  required  at  the  upstream 
boundary  to  keep  the  discharge  constant  at  various  reser¬ 
voir  water  levels  is  then  determined. 

A  dimensionless  curve,  which  is  a  plot  between 
H  /H  and  d/H  for  a  constant  discharge  is  developed 

i  L|  III  a  a  III  u  X 

and  presented  on  Fig.  3.7.  In  this  plot  Hmax  is  the  initial 

head  at  the  upstream  boundary  corresponding  to  the  condition 

when  there  is  no  water  in  the  reservoir.  H^  represents 

the  head  required  at  the  upstream  boundary  to  give  the 

same  discharge  at  various  reservoir  water  levels  d.  The 

line  of  constant  discharge  shows  that  this  relationship  is 

also  highly  non-linear.  It  may  be  seen  further  that  the 

increase  in  head  required  at  low  reservoir  water  levels  is 

not  very  significant,  whereas  at  the  higher  water  levels 

the  increase  in  head  required  is  quite  appreciable,  e.g., 

for  ratios  of  d/H  from  0.8  to  1.0,  the  increase  in  head 

max 

required  ranges  from  20  to  40%  (Fig.  3.7). 

This  curve  may  be  used  for  determining  the  head  re¬ 
quired  at  various  reservoir  water  levels  for  a  typical 
value  of  q/kll  .  From  this  curve  the  head  required  at  the 
upstream  boundary  for  various  reservoir  water  levels  is 
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FIG  3  7  DIMENSIONLESS  PLOT  OF  HEAD  REQUIRED  FOR 
DIFFERENT  RESERVOIR  LEVELS 
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determined  and  then  the  free  surface  and  potentials  are 
obtained  for  the  co rres pondi ng  reservoir  levels  ranging 
from  400  to  1000  in  increments  of  200.  As  the  increase 
in  head  at  200  elevation  of  the  reservoir  water  level  is 
very  small,  this  case  has  not  been  studied. 

These  plots  of  free  surface  and  eq ui po tenti a  1 s  are 
shown  on  Fig.  3.8  through  3.11. 
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3 . 4  Effect  of  Changes  in  Reservoir  Level  on  the  Stability 

Ignoring  the  Terzaghi  and  Skempton  softening  mechanism, 
almost  all  the  studies  dealing  with  the  long  term  stability 
of  natural  slopes  and  embankment  based  on  slip  circle 
analysis  show  that  the  factor  of  safety  increases  with  the 
raising  of  water  level  in  the  reservoir.  The  horizontal 
component  of  the  hydrostatic  pressure  obtained  as  a  result 
of  the  filling  of  reservoir  stabilizes  the  slopes  and  as 
such  with  the  raising  of  water  level  in  the  reservoir  the 
factor  of  safety  increases.  This  observation  is  however 
based  on  the  assumption  that  the  raising  of  the  water  level 
in  the  reservoir  does  not  influence  the  head  at  the  up¬ 
stream  boundary. 

It  is  already  established  in  Section  3.3  that  the 
impounding  of  the  reservoir  does  influence  the  head  at  the 
upstream  boundary  and  therefore,  it  is  not  necessary  that 
the  factor  of  safety  should  improve  at  higher  reservoir 
water  levels.  The  increase  in  head  at  the  upstream  boundary 
corresponding  to  the  respective  reservoir  water  levels  re¬ 
sults  in  an  increase  of  pore  pressure  in  the  soil  or  rock 
mass,  which  in  turn,  reduces  the  stability. 

To  assess  the  implications  of  this  hypothesis  on  the 

stability  of  slopes,  a  typical  30°  slope  with  an  impermeable 

base,  consisting  of  homogeneous  and  isotropic  material  hav- 

2 

ing  strength  parameters  of  <j> 1  =  40°,  c1  =  5 0 1 / m  ,  y 
2.2t/m^  is  analyzed. 

As  shown  on  Fig.  3.12  it  may  be  seen  that  the  factor 
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FIG.  3.12  VARIATION  OF  FACTOR  OF  SAFETY  WITH  RESERVOIR  LEVELS 
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of  safety  obtained  by  using  the  Bishop  simplified  method 
(Bishop  (1955))  increases  with  the  raising  of  water  level 
in  the  reservoir,  provided  the  head  at  the  upstream  bound¬ 
ary  is  kept  constant.  If  the  head  at  the  upstream  boundary 
is  adjusted  in  accordance  with  the  reservoir  level  in  order 
to  keep  discharge  constant  and  if  the  pore  pressure  data 
thus  obtained  with  these  boundary  conditions  are  used  in 
the  stability  analysis,  it  may  be  seen  that  the  factor  of 
safety  increases  up  to  a  reservoir  level  650  and  then  starts 
decreasing  with  further  raising  of  the  reservoir  level 
(Fig.  3.12). 

It  is  also  evident  that  this  effect  is  more  signifi¬ 
cant  at  higher  reservoir  levels,  which  is  in  contrast  to 
most  of  the  existing  conventional  studies.  A  reduction  of 
about  9%  in  the  factor  of  safety  is  witnessed  when  the  re¬ 
servoir  level  is  at  1000  elevation.  At  lower  reservoir 
levels  the  reduction  is  however  not  very  significant  but  may 
be  significant  if  the  slope  is  already  in  a  limiting  state. 

To  illustrate  the  importance  of  this  reduction  in 
factor  of  safety  more  explicitly  a  case  of  reservoir  level 
elevation  at  1000.0  is  selected  and  using  different  strength 
parameters  a  comparison  of  factor  of  safety  obtained  with  the 
two  hydraulic  boundary  conditions  is  made.  Three  values  of 
angle  of  shearing  resistance  (^')  namely  40°,  30°  and  20° 
are  selected  and  for  each  value  of  <J> '  ,  three  values  of  co¬ 
hesion  ( c 1 )  ,  i.e.,  50t/m2,  25t/m2  and  5t/m2  are  tried.  A 
relationship  between  factor  of  safety  and  angle  of  shearing 


' 


62 


resistance  is  obtained  for  different  values  of  c1  and 
curves  are  developed  using  both  the  constant  head  and  the 
constant  discharge  methods  (Fig.  3.13  through  3.15). 

As  evident  from  Fig.  3.13  a  typical  slope  of  soil  or 

P 

rock  mass  having  strength  parameters  c'  =  50t/m  and 
<t> 1  =  29°  would  be  considered  as  marginally  stable  if  the 
factor  of  safety  is  obtained  by  using  the  constant  head 
condition,  whereas  on  the  basis  of  the  constant  discharge 
concept,  the  slope  would  fail. 

Similarly,  referring  to  Fig.  3.14  and  Fig.  3.15  slopes 

p 

having  strength  parameters  as  c'  =  25t/m  ,  <j> 1  =  31°  and 
c‘  =  5t/m  ,  <p‘  =  33  1/2°  respectively  may  be  treated  as 
marginally  safe  under  constant  head  conditions  but  are  un¬ 
stable  when  the  constant  discharge  conditions  are  invoked. 

This  obviously  shows  that  even  if  the  reduction  in  factor 
of  safety  is  small,  in  the  case  of  slopes  of  marginal  sta¬ 
bility,  the  effect  of  the  constant  discharge  concept  might 
be  found  very  crucial. 

Studies  have  already  shown  that  when  dealing  with  the 
long  term  stability  of  slopes,  especially  of  fissured  clay 
or  rocks,  a  reduction  in  cohesion  takes  place  due  to  a  soften¬ 
ing  mechanism,  which  may  be  aggravated  by  the  submergence  of 
the  slopes  (Skempton  (1964);  DeLory  (1957)). 

If  this  phenomenon  is  taken  into  account  the  factor 
of  safety  for  a  30°  slope  having  values  of  c‘  and  <p 1  as 
50t/m2  and  33  1/2°  respectively  obtained  as  1.25.  Assuming 

p 

the  value  of  c'  is  reduced  to  5t/m  in  due  course  of  time, 
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FIG.  3.13  RELATIONSHIP  BETWEEN  FACTOR  OF  SAFETY  AND  ANGLE  OF  SHEARING  RESISTANCE 
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FIG.  3.14  RELATIONSHIP  BETWEEN  FACTOR  OF  SAFETY  AND  ANGLE  OF  SHEARING  RESISTANCE 
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FIG.  3.15  RELATIONSHIP  BETWEEN  FACTOR  OF  SAFETY  AND  ANGLE  OF  SHEARING  RESISTANCE 
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it  may  be  seen  that  the  slope  is  found  to  be  marginally 
safe  with  the  constant  head  condition  but  if  the  stability 
analysis  is  based  on  the  constant  discharge  concept,  the 
same  slope  is  unstable  (Fig.  3.13  through  3.15}. 

However,  other  things  being  equal,  i.e.,  the  c'  and 
<j> 1  values  remain  constant,  the  impounding  results  in  an 
increase  in  factor  of  safety  when  slip  surface  is  an  arc 
of  a  ci rcl e,  regard! ess  of  the  discharge  conditions  (Fig. 

3.12).  This  increase  is  less  in  the  case  of  constant  discharge 
than  that  of  constant  head  conditions. 

It  is  also  evident  that  when  there  is  no  water  in 
the  reservoir  a  factor  of  safety  of  1.358  is  obtained  with 
c'  =  50t/m^  and  <j> 1  =  40°.  As  the  filling  proceeds,  the 
factor  of  safety  increases.  At  the  1000  reservoir  elevation 
if  it  is  assumed  that  the  cohesion  is  reduced  to  5t/m  , 
the  factor  of  safety  obtained  by  the  constant  head  condi¬ 
tions  is  1.411  which  is  still  higher  than  the  initial  factor 
of  safety  correspondi ng  to  no  water  in  the  reservoir.  On 
the  contrary,  if  the  constant  discharge  concept  is  used  for  de¬ 
termining  the  factor  of  safety  at  the  1000  reservoir  elevation 

with  c'  =  5t/m2  and  <J> '  =  40°,  a  factor  of  safety  of  1.246 
is  obtained  which  is  lower  than  the  initial  factor  of  safety 
and  thus  shows  that  this  condition  may  trigger  a  slide. 

It  is  also  interesting  to  note  that  in  this  typical 
example  of  a  slope,  extreme  conditions  are  not  considered; 
that  is  to  say  the  initial  slope  of  the  ground  water  table 
before  impounding  is  about  1  in  5,  the  upstream  head  is 
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taken  at  an  elevation  of  1000  and  its  location  is  chosen 
at  a  distance  four  times  the  head  inside  the  slope.  There 
could  be  more  unfavourable  conditions  which  one  can  en¬ 
counter  in  the  field.  As  an  example  it  may  be  quoted  that 
depending  upon  the  geo  hydro! ogi cal  conditions  of  the  region 
and  the  local  topography,  the  initial  slope  of  the  water 
table  may  be  still  steeper.  This  would  simply  mean  that 
the  amount  of  outflow  from  the  banks  would  be  more,  which 
in  turn,  results  in  a  further  increase  in  head  required 
at  the  upstream  boundary  for  various  reservoir  levels.  In 
that  case  a  further  increase  in  pore  pressure  due  to  these 
hydraulic  boundary  conditions  would  make  the  reduction  of 
factor  of  safety  much  more  significant. 

It  may  further  be  noted  that  if  non-homogeneity  of 
the  soil  or  rock  mass  is  considered  in  conjunction  with  the 
constant  discharge  conditions,  the  initial  head  at  the 
u/s  boundary  may  further  rise,  resulting  in  higher  pore 
pressures  and  thus  a  further  reduction  in  the  factor  of 
safety.  No  detailed  analysis  for  non-homogeneous  material 
has  been  done,  yet  qualitatively  the  effect  of  non -homogenei ty 
can  be  understood  from  a  hypothetical  model  shown  on  Fig. 

3.16. 

As  shown  on  this  figure  if  the  flow  is  occurring 
from  low  to  high  permeability  material,  since  the  discharge 
and  the  total  head  at  the  interface  is  the  same,  the  free 
surface  line  in  the  less  permeable  material  will  be  more 
steeper.  As  soon  as  it  enters  the  more  permeable  material. 


FREE  .SURFACE  LINE  FOR  NON-HOMOGENEOUS  MATERIAL 
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the  free  surface  line  becomes  flatter  and  if  the  difference 
in  the  two  permeabilities  is  large, the  flow  in  the  lower 
layer  near  the  interface  will  almost  be  horizontal.  There¬ 
fore,  it  is  evident  that  under  these  conditions  at  the 
chosen  location  of  the  upstream  boundary,  the  head  will 
further  rise  due  to  this  non-homogeneity,  producing  thereby 
an  additional  adverse  effect  on  the  stability. 
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3 . 5  Influence  of  Infiltration  on  the  Stability 

To  evaluate  the  influence  of  infiltration  on  the 
stability  of  slopes,  a  constant  infiltration  of  5  cm/ day 
over  the  entire  free  surface  is  taken.  A  dimensionless 
infiltration  factor  (Q)  is  introduced  which  is  defined 
as  the  volume  of  water  infiltrating  per  unit  horizontal 
area  per  unit  time  divided  by  the  coefficient  of  permeabi¬ 
lity  of  the  soil  or  rock  mass  (i.e.,  Q  =  -  q/k).  The  co¬ 
efficient  of  permeability  is  taken  as  100  cm/day,  i.e., 

_  3 

1.1  x  10  cms/sec.  The  computer  program  is  used  to  deter¬ 
mine  the  free  surface  and  potentials  for  the  case  when  the 
reservoir  level  is  at  800  (Fig.  3.18).  The  change  in  the 
free  surface  due  to  this  infiltration  is  shown  in  Fig  3.17 
and  Fig.  3.18.  As  evident  this  results  in  the  raising  of 
the  free  surface  line.  The  pore  pressure  data  obtained 
from  these  flow  conditions  is  then  used  in  the  stability 
analysis. 

The  case  of  upstream  head  at  1200  elevation  and  re¬ 
servoir  level  at  800  with  c1  =  50t/m3,  <J> 1  =  40°  and  no 
infiltration  over  the  free  surface  gives  a  factor  of  safety 
of  1.471.  When  a  constant  infiltration  of  5  cm/day  is  in¬ 
troduced  and  assuming  the  coefficient  of  permeability  as 
1.1  x  10~3  cm/sec,  the  factor  of  safety  obtained  is  1.439. 
This  obviously  shows  that  the  influence  of  direct  infiltra¬ 
tion  on  the  increase  in  pore  pressures  and  thus  on  the 
stability  of  slopes  is  not  very  significant.  This  is  only 
true  when  the  effect  of  infiltration  is  considered  in  terms 
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FIG.  3.17  DIAGRAM  SHOWING  INFLUENCE  OF  INFILTRATION  ON  THE  FREE  SURFACE 
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FIG.  3.18  FREE  SURFACE  AND  EQUIPOTENTIALS  FOR  TYPICAL  SLOPE  (30  ) 
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of  raising  of  free  water  table  only  which  in  turn  results 
in  the  increase  of  pore  pressures  in  the  stability  analysis. 

No  effect  of  infiltration  on  the  filling  of  cracks  result¬ 
ing  in  the  joint  water  pressure  or  softening  etc.  is  taken 
into  consi deration. 

Infiltration  over  the  entire  free  surface  as  considered 
here,  is  a  direct  consequence  of  precipitation  which  depends 
not  only  on  the  geographic  location  but  on  the  permeability 
of  the  material  as  well.  It  may  be  seen  that  with  even  more 
favourable  values  of  run  off  factor  and  permeability  an 
infiltration  of  5  cm/day  over  the  entire  free  surface  may 
only  be  obtained  from  very  high  values  of  precipitation. 

It  may  therefore  be  concluded  that  barring  very  ex¬ 
treme  conditions  of  precipitation,  direct  infiltration  on 
the  free  surface  will  result  in  the  raising  of  the  free 
surface  line.  This  rise  in  the  free  ground  water  table 
due  to  infiltration  does  not  induce  a  very  significant  in¬ 
increase  in  the  pore  pressures.  Therefore,  ignoring  the 
effect  of  filling  of  cracks  by  water  and  softening  of  material 
above  the  free  ground  water  table,  this  increase  in  pore 
pressures  hardly  produces  any  significant  effect  on  the 
stability  of  slopes. 
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CHAPTER  IV 


DISCUSSION  OF  RESULTS  AND  CONCLUSIONS 

The  present  work  has  been  carried  out  with  the  object 
of  evaluating  the  effect  of  changes  in  reservoir  level  on 
the  stability  of  natural  slopes.  Although  a  detailed 
account  of  the  results  obtained  is  given  in  the  preceding 
chapters,  a  brief  summary  of  the  results  and  conclusions 
drawn  is  given  in  the  following  paragraphs. 

(1)  A  Fortran  IV  computer  program  using  finite  differ¬ 
ence  techniques  has  been  developed  which  can  be  used  for 
determining  the  steady  state  free  surface  potential  flow 
through  porous  media.  The  program  is  quite  versatile  and 

can  be  used  for  any  geometric  configuration  and  boundary 
conditions.  As  a  check,  the  numerical  solution  obtained 
by  using  this  program  is  compared  with  Numerov's  analytical 
solution  of  free  surface  seepage  through  a  typical  dam 
(Fig.  1.2).  A  close  agreement  between  the  two  results  is 
obtained.  A  value  of  1.6  for  the  over-relaxation  factor 
ensured  the  quickest  convergence.  This  program  can  also 
be  used  for  solving  free  surface  problems  with  constant 
infiltration.  This  has  been  used  in  determining  the  free 
surface  and  potentials  when  a  typical  constant  infiltration 
is  introduced  over  the  entire  free  surface. 

(2)  A  study  for  the  outflow  from  the  banks  at  various 
reservoir  levels  shows  that  the  relationship  between  the 
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reservoir  level  and  the  discharge  is  highly  non-linear. 

At  lower  reservoir  levels  (up  to  400  elevation),  the  re¬ 
duction  in  discharge  is  about  6  % ,  but  at  higher  levels, 
e.g.,  at  800  elevation,  the  reduction  in  discharge  is  almost 
50%  (Fig.  3.6).  This  is  however  based  on  the  assumption 
that  the  head  at  the  upstream  boundary  remains  constant. 

(3)  A  hypothetical  concept  of  constant  discharge 
is  postulated  which  results  in  a  raising  of  head  at  the 
upstream  boundary  with  the  raising  of  reservoir  level. 

The  occurrence  and  movement  of  ground  water  is  an 
important  phase  of  the  hydrological  cycle  which  always  main¬ 
tains  a  balance  "of  inflow  and  outflow  in  the  ground  water 
basin.  Any  change  in  the  exit  conditions  at  the  face  of 
the  slope  which  may  correspond  to  the  filling  of  reservoir 
need  not  result  in  the  reduction  of  discharge  from  the  bank 
especially  in  the  areas  of  high  relief.  This  obviously 
means  that  the  assumption  of  constant  head  at  the  upstream 
boundary  which  is  generally  considered  to  be  uninfluenced 
by  the  geometric  configuration  and/or  the  exit  conditions 
of  the  slope  is  at  variance  with  reality.  It  is  therefore 
concluded  that  in  defining  the  flow  domain  under  steady 
state  seepage  conditions  the  concept  of  constant  discharge 
which  results  in  an  increase  in  head  at  the  upstream  boun¬ 
dary  with  the  raising  of  reservoir  level  is  more  realistic 
for  field  cases. 

The  transient  conditions  would  always  be  on  the  safe 
side  as  far  as  the  stability  of  slopes  is  concerned. 


. 
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It  is  observed  that  the  relationship  between  the 
head  required  at  the  upstream  boundary  to  keep  the  discharge 
constant  and  the  reservoir  level  is  also  non-linear  (Fig. 
3.7).  At  lower  reservoir  levels  the  increase  is  not  very 
significant  but  at  higher  levels  the  increase  in  head  re¬ 
quired  is  quite  appreciable,  e.g.,  at  1000  reservoir  eleva¬ 
tion,  the  increase  is  40%.  Quantitatively  this  increase  in 
head  at  the  upstream  boundary  due  to  raising  of  reservoir 
level  results  in  a  raising  of  the  free  surface  line  and 
thereby  an  increase  in  the  pore  pressures  in  the  stability 
analyses . 

(4)  The  stability  of  a  typical  30°  slope  consisting 
of  a  homogeneous,  isotropic  material  and  having  an  imperme¬ 
able  base  is  analyzed  at  various  reservoir  levels  ranging 
from  0  to  1000  elevation  in  increments  of  200. 

The  effect  of  impounding  of  the  reservoir  on  the 
stability  of  slopes  is  then  compared  by  using  the  constant 
head  and  the  constant  discharge  methods.  The  material  pro¬ 
perties  of  c‘  =  50t/m^,  <j> 1  =  40°  and  y^  =  2.2t/m^  are  used 
in  the  Bishop  simplified  method  for  analyzing  the  stability 
of  this  typical  slope.  A  minimum  factor  of  safety  is 
obtained  in  each  case  corresponding  to  various  reservoir 
levels  and  hydraulic  boundary  conditions  at  the  upstream 
boundary.  They  are  given  in  Table  4.1  and  Fig.  3.12. 

It  is  observed  that  the  factor  of  safety  increases 
with  the  raising  of  reservoir  level  under  constant  head 
conditions.  The  factor  of  safety  increases  up  to  a  reservoir 
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TABLE  4.1 

VARIATION  OF  FACTOR  OF  SAFETY  WITH  RESERVOIR 


LEVEL  AND 

UPSTREAM 

BOUNDARY 

HEAD 

c 1  = 

50t/m2 

4'  =  40 

o 

Head  At 

The  Upstream 
Boundary 

Min  . 

Factor 

of  Safety 

for  Re  servo ir 

Level 

0 

200 

400 

600 

800 

1000 

1000 

1  .358 

1  .443 

1  .484 

1  .518 

1  .540 

1  .556 

1055 

1  .453 

1100 

1  .483 

1200 

1  .471 

1  400 

1  .424 

1200 

with 

Q  =  0.05 

1  .439 

Q  =  Infiltration  Factor 
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elevation  at  650  and  then  starts  decreasing  under  constant 
discharge  conditions.  A  decrease  in  the  factor  of  safety 
of  about  9%  is  witnessed  at  reservoir  elevation  1000  if  the 
constant  discharge  conditions  are  invoked.  However,  the 
factor  of  safety  still  remains  higher  than  the  initial 
value  of  1.358  corresponding  to  no  water  in  the  reservoir. 

This  shows  that  in  this  particular  case  the  filling  of  re¬ 
servoir  produces  a  stabilizing  effect  even  under  the  con¬ 
stant  discharge  conditions. 

It  may,  however,  be  noted  that  in  this  particular 
problem  extreme  conditions  are  not  considered.  The  amount 
of  initial  outflow  depending  on  the  initial  slope  of  the 
ground  water  table  is  such  that  at  reservoir  elevation  1000, 
a  4 0 %  increase  in  the  head  at  the  upstream  boundary  is  re¬ 
quired.  Depending  upon  the  geohydrol ogi cal  conditions  and 
the  local  topography,  there  could  be  field  cases  when  the 
initial  slope  of  the  ground  water  table  is  steeper.  In 
that  case  if  the  constant  discharge  concept  is  used,  the 
head  required  at  the  upstream  boundary  will  follow  a  differ¬ 
ent  relationship  and  obviously  higher  values  will  be  obtained. 
This  means  the  rise  in  the  ground  water  table  will  be  more 
which  will  induce  higher  pore  pressures.  Under  these 
conditions  it  is  possible  that  at  higher  reservoir  levels 
the  factor  of  safety  may  become  less  than  the  initial  value, 
and  thereby  may  initiate  instability. 

It  is  observed  that  in  the  cases  of  marginal  stability 
the  application  of  the  constant  discharge  concept  results 
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in  instability  (Fig.  3.13  through  Fig.  3.15). 

(5)  To  investigate  further  effects  of  the  constant 
discharge  concept,  a  case  of  the  reservoir  level  at  1000 
elevation  and  upstream  head  at  1400  elevation  is  studied 
in  detail.  Different  combinations  of  c'  and  <j> 1  [e.g., 
c1  =  50t/m^,  25t/m^,  5t/m^  and  <J> 1  =  40°,  30°  20°)  have  been 
tried  and  sets  of  curves  showing  the  relationship  between 
the  factor  of  safety  and  the  angle  of  shearing  resistance 
are  obtained  (Fig.  3.13  through  Fig.  3.15).  The  minimum 
factor  of  safety  obtained  by  using  the  constant  head  and 
the  constant  discharge  conditions  for  different  combinations 
of  c1  and  <f> 1  are  presented  in  Table  4.2. 

2 

It  is  observed  that  a  typical  slope  with  c1  =  50t/m 
and  (p '  =  40°  has  a  factor  of  safety  of  1.358  before  the 
filling  of  the  reservoir.  Under  the  constant  head  condi¬ 
tions  the  factor  of  safety  increases  with  the  raising  of 
the  reservoir  level  even  if  the  value  of  c‘  is  reduced  to 

p 

5t/m  .  A  factor  of  safety  of  1.411  is  obtained  for  c  = 

o 

5t/m  and  <j> 1  =  40°  under  constant  head  conditions,  which 
is  still  higher  than  the  initial  value  of  1.358  (Table  4.2). 
If  the  constant  discharge  conditions  are  considered,  it 
is  observed  that  a  factor  of  safety  of  1.246  is  obtained 
for  c'  =  5t/m2  and  <j> '  =  40°.  This  obviously  shows  that  the 
factor  of  safety  has  dropped  down  from  its  initial  value 
and  thus  may  result  in  triggering  a  slide. 

An  overall  conclusion  can  therefore  be  drawn  that 
depending  upon  the  geohydrol ogi cal  conditions,  local  topo- 
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TABLE  4.2 

COMPARISON  OF  FACTOR  OF  SAFETY 


Reservoi r 

Elevation 

=  1000 

( Degrees ) 

c1  =  50t/m^ 

c1  =  25t/m^ 

c'  = 

5  t/m^ 

F.S.(l) 

F.S. (2) 

F.S.(l) 

F.S . (2) 

F.S.(l) 

F.S. (2) 

40 

1  .556 

1  .424 

1  .491 

1  .349 

1.411 

1  .246 

30 

1  .  10 

1  .022 

1  .047 

0.953 

0.977 

0.87 

20 

0.727 

0.687 

0.682 

0  .627 

0.623 

0.560 

F.S.(l)  = 

Factor  of 

safety  with  constant  head. 

F.S.(2)  =  Factor  of  safety  with  constant  discharge. 

Upstream  Head  =  1000  (under  constant  head  conditions). 
Upstream  Head  =  1400  (under  constant  discharge  conditions). 
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graphy  and  the  material  properties  a  stable  slope  corres¬ 
ponding  to  no  water  in  the  reservoir,  which  is  usually  con¬ 
sidered  the  worst  case  may  become  unstable  at  higher  re¬ 
servoir  elevations  if  the  constant  discharge  conditions  are 
invoked.  Therefore,  it  is  important,  especially  when  work¬ 
ing  in  the  areas  of  high  relief,  to  have  a  good  knowledge 
of  the  geohydrol ogi cal  conditions  of  the  area,  so  that  more 
realistic  hydraulic  boundary  conditions  may  be  adopted  in 
defining  the  flow  domain. 

(6)  The  stability  analyses  are  done  only  for  homoge¬ 
neous,  isotropic  material.  Although  no  quantitative  ex¬ 
planation  is  given,  it  is  of  interest  to  note  that  the 
effect  of  non-homogeneity  in  addition  to  the  constant  dis¬ 
charge  conditions  produces  an  even  more  adverse  effect  on 
the  stability  of  slopes  (Fig.  3.16). 

(7)  The  effect  of  direct  infiltration  on  the  free 
surface  and  then  on  the  pore  pressures,  which  are  used  in 
the  stability  analysis  is  evaluated.  A  case  of  reservoir 
level  at  800  elevation  and  upstream  boundary  head  at  1200 

elevation  with  c'  =  50t/m  and  <j> '  =  40°  is  investigated. 

_  3 

The  coefficient  of  permeability  is  taken  as  1.1  x  10  cm/sec. 
It  is  found  that  in  this  typical  case  the  effect  of  infil¬ 
tration,  considered  only  in  terms  of  an  increase  of  pore 
pressure  affects  only  the  third  significant  figure  in  the 
factor  of  safety. 
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CHAPTER  V 


POSSIBLE  EXPLANATION  OF  VAJONT  ROCK  SLIDE 

On  the  basis  of  the  constant  discharge  concept,  it 
is  now  possible  to  give  a  qualitative  explanation  of  the 
sliding  movements  at  Vajont  which  occurred  during  the 
first  stage  of  filling.  Muller  (1964)  has  mentioned  that 
a  few  springs  were  observed  at  the  toe  of  the  left  slope  in 
the  Vajont  gorge  before  filling  the  reservoir.  The  piezo¬ 
metric  observations  in  two  piezometers  were  started  in 
August,  1961  and  in  the  third  piezometer  in  October,  1961 
(Muller,  1964). 

At  piezometer  (PZ^)  located  at  the  farthest  point 
from  the  face  of  the  slope  the  observations  which  were 
started  in  October,  1961  showed  a  water  level  close  to  700 
m.  The  other  two  piezometers  (PZ-j  and  PZ^)  which  were 
located  close  to  the  valley  recorded  water  levels  of  about 
630  m.  The  reservoir  level  was  fluctuating  around  590  m. 
at  that  time. 

These  observations  establish  the  fact  that  at  the  out 
set  the  mountain  water  table  was  sloping  towards  the  valley 
and  there  was  an  outflow  from  the  slopes  into  the  valley. 

As  no  piezometric  observations  during  the  first  stage  of 
filling  are  available,  the  slope  of  the  mountain  water  tabl 
was  not  known  at  that  time.  During  the  first  filling  when 
the  reservoir  level  reached  elevation  650  m.  accelerated 
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movements  started.  The  reservoir  level  was  then  lowered 
and  at  elevation  600  m.  the  movement  practically  stopped. 

Until  November,  1961  the  lake  water  level  was  kept  constant 
and  there  was  practically  no  displacement.  From  the  beginn¬ 
ing  of  October,  1961  to  the  beginning  of  February,  196?., 
the  water  level  was  raised  from  590  m.  to  650  m.  and  the 
movement  remained  small.  Rocks  started  to  move  again  as 
soon  as  the  water  level  rose  above  the  level  650  m .  These 
movements  continued  until  the  beginning  of  1963,  when  a 
short  slow  down  occurred  as  the  reservoir  level  was  again 
dropped  from  700  m.  to  about  650  m.  During  the  third  fill¬ 
ing  from  the  beginning  of  April  to  the  end  of  May,  1963, 
the  lake  level  rose  from  650  m.  to  696  m.  and  the  displace¬ 
ments  were  again  found  negligible.  When  the  water  rose 
again  from  696  m.  to  710  m.  rock  displacements  became  alarm¬ 
ing  and  resulted  in  the  final  phase  of  the  failure. 

Muller  (1968)  has  tried  to  explain  this  first  sliding 
movement  by  considering  an  artesian  thrust  of  about  100000  t/m 
at  the  base  of  the  sliding  mass.  Although  nothing  was  known 
about  the  permeability  of  the  individual  strata,  Muller 
(1968)  has  assumed  that  the  permeability  of  the  strata  ly¬ 
ing  under  the  slip  surface  (Dogger)  was  more  than  the  strata 
above,  which  resulted  in  this  hydrodynamic  pressure  with 
a  vertical  component,  correspondi ng  to  a  maximum  height  of 
water  column  of  150  m.  He  further  explains  that  during  the 
first  movement  the  friction  descreased  from  its  peak  value 
to  the  smaller  residual  values.  Yet  the  reduced  frictional 
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resistance  increased  again  when  a  state  of  rest  sets  in, 
the  contact  material  consolidated  and  when  new  movements 
started  the  old  peak,  strength  value  (or  value  near  to  it) 
was  prevailing. 

If  Muller's  assumptions  that  the  material  underlying 
the  slip  surface  was  more  permeable  and  that  the  artesian 
thrust  caused  the  first  sliding  movement  at  the  reservoir 
level  635  m.  is  accepted,  during  the  second  filling  the 
same  phenomenon  should  have  occurred  at  the  reservoir  level 
635  m.  But  no  movements  were  observed  until  the  reservoir 
elevation  reached  650  m.  His  explanation  of  regaining 
strength  when  a  state  of  rest  sets  in  is  at  variance  with 
real i ty . 

Kenney  (1967)  is  of  the  opinion  that  the  reservoir 
slope  must  have  been  in  a  very  delicate  state  of  stability 
for  a  long  period  of  time.  Assuming  the  material  as  purely 
frictional,  he  has  found  that  the  factor  of  safety  decreases 
with  the  filling  of  the  reservoir  up  to  900  m.  elevation. 

At  700  m.  elevation  the  decrease  is  5-10%  while  at  635  m. 
elevation  the  decrease  is  only  4%.  Although  this  decrease 
v/ as  small  it  caused  the  rock  to  move.  He  also  assumes  that 
since  the  coefficient  of  sliding  resistance  was  less  than 
0.4  ( correspondi  ng  to  <j> '  =  22°)  when  the  reservoir  slope 
failed,  a  large  portion  of  the  sliding  surface  passed 
through  the  clay  material  rather  than  through  limestone. 

A  critique  of  Kenny's  explanation  is  given  as  follows: 

The  assumption  that  the  rock  material  was  purely  frictional 
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is  rather  debateable.  Therefore  it  is  not  necessary  that 
the  coefficient  of  sliding  resistance  may  decrease  with 
the  filling  of  reservoir  as  hypothetically  postulated  by 
him. 

Moreover,  even  if  it  is  assumed  that  the  slope  material 
is  purely  frictional  and  that  the  sliding  surface  passed 
through  the  clay  layers  which  had  residual  strength,  the 
first  filling  of  reservoir  up  to  635  m.  elevation  resulted 
in  the  reduction  of  the  coefficient  of  sliding  resistance 
by  about  4%  and  caused  the  rock  slope  to  move.  But  during 
the  second  filling  no  movements  were  observed  until  the 
reservoir  level  reached  the  650  m.  elevation.  Therefore, 
it  is  argued  that  if  the  material  is  considered  frictional 
and  was  at  residual  state,  the  second  filling  of  the  reser¬ 
voir  at  elevation  635  m.  would  have  again  given  the  previous 
value  of  the  mobilized  coefficient  of  the  sliding  resistance 
(i.e.,  Mo  reduction  in  the  initial  coefficient  of  sliding 
resistance  at  635  m.  reservoir  elevation).  The  movements 
should  have  then  started  at  635  m.  elevation  again  and  not 
at  650  m.  elevation  during  the  second  stage  of  filling. 

On  the  basis  of  the  low  value  of  angle  of  shearing 
resistance  obtained  at  failure  Kenney  (1967)  concluded  that 
the  slip  surface  passed  through  the  clayey  material  which 
was  at  residual  state.  Muller  (1968)  contradicted  this 
possibility  and  in  support  has  given  results  of  detailed 
geological  investigations  showing  thin  layers  of  clayey 
material  only  very  sparsely  present  between  the  bedding 
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planes.  Hence  the  possibility  of  the  slip  surface  passing 
through  these  thin  layers  of  clayey  material  is  remote. 

Jaeger  (1969)  has  explained  this  first  time  slide 
by  assuming  low  values  of  angle  of  shearing  resistance  along 
the  upper  strata.  He  assumes  that  during  this  period  the 
factor  of  safety  was  very  close  to  unity  and  when  the  water 
receded  it  became  more  than  one.  The  disappearance  of  any 
movement  during  the  second  stage  of  filling  up  to  the  former 
maximum  level  of  658  m.  is  attributed  to  the  remarkable 
process  of  "auto-stabi 1 i zation" . 

It  can  be  argued  that  no  field  evidence  is  available 
to  prove  the  assumption  that  the  angle  of  shearing  resist¬ 
ance  of  the  upper  strata  was  lower  than  the  lower  strata. 

On  the  contrary  it  seems  more  reasonable  to  assume  lower 
values  of  friction  in  the  lower  strata  which  was  submerged. 
Moreover  the  process  of  auto -s tabi 1 i za ti on  as  given  by  him 
does  not  explain  fully  the  reason  as  to  why  the  displacements 
did  not  start  again  at  the  reservoir  elevation  635  m.  dur¬ 
ing  the  second  stage  of  filling.  If  the  material  was  at 
this  residual  state,  the  displacements  should  have  started 
again  at  reservoir  elevation  635  m.  as  no  regain  of  strength 
is  possible  even  by  this  process  of  auto-stabilization. 

To  summarize  it  may  be  stated  that  although  many  ex¬ 
cellent  references  are  available  which  describe  the  possible 
causes  of  the  Vajont  rock  slide,  none  of  them  is  free  from 
unrealistic  assumptions  and  unresolved  questions. 

In  spite  of  the  scanty  information,  the  rock  displace- 
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merits  in  the  Vajont  valley  corresponding  to  the  various 
reservoir  levels  can  be  explained  qualitatively  on  the 
basis  of  the  constant  discharge  concept. 

It  may  be  argued  that  before  the  beginning  of  the 
first  filling,  the  mountain  water  table  was  steep.  The 
initial  factor  of  safety  correspond!*  ng  to  the  case  of  no 
impounding  must  have  been  more  than  unity.  As  the  filling 
proceeded  the  water  table  must  have  risen  upward  in  order 
to  keep  the  discharge  constant.  Due  to  this  rise  in  mountain 
water  table  pore  pressure  increased  which  in  turn  reduced 
the  factor  of  safety  very  close  to  unity  and  the  rock  slopes 
started  to  move.  During  the  period  of  lowering  and  then 
of  constant  reservoir  level,  the  movement  stopped  as  the 
water  table  came  down  and  the  factor  of  safety  become  more 
than  one  again.  When  the  second  filling  started  the  rock 
slope  did  not  move  until  the  former  reservoir  elevation  of 
650  m.  was  crossed.  This  is  quite  evident  from  the  fact 
that  during  the  second  filling  the  water  table  started  ris¬ 
ing  again  and  as  soon  as  it  reached  to  the  previous  posi¬ 
tion  corres pondi ng  to  the  reservoir  level  650,  the  factor 
of  safety  became  very  close  to  unity  and  the  rock  slope 
started  to  move  again. 

It  is  of  interest  to  note  that  the  rock  material  which 
was  submerged  became  more  permeable  and  during  the  last 
stages  of  filling  exhibited  fully  drained  conditions.  It 
is  expected  that  the  rock  material  in  the  beds  stratified 
horizontally  became  buoyant  due  to  submergence  and  as  a 
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result  of  this  stress  release,  the  joints  and  fissures 
opened  up.  Hence  the  permeability  of  the  rock  material  in 
these  horizontal  layers  increased,  w hereas  the  permeability 
of  the  inclined  layers  of  Malm  and  lower  cretaceous  forma¬ 
tions  which  were  not  submerged  remained  the  same.  As  the 
amount  of  water  flowing  out  of  the  ground  water  basin  is 
constant  and  the  head  loss  in  the  less  permeable  material 
is  more,  it  becomes  necessary  that  the  water  table  must 
have  a  steeper  gradient  in  these  inclined  layers  of  com¬ 
paratively  low  permeabi 1 i ty .  Hen  ce  it  is  con  cl uded  that 
due  to  this  non -homogenei ty  the  slope  of  the  mountain  water 
table  in  the  rock  beds  stratified  obliquely  and  horizontally 
was  not  the  same. 

With  the  filling  of  the  reservoir,  since  the  exit 
conditions  are  disturbed  and  the  excess  of  water  in  the 
ground  water  basin  has  to  get  out,  the  effect  of  raising 
the  reservoir  on  the  water  table  in  the  less  permeable 
material  will  be  more  adverse.  The  water  table  in  the  more 
permeable  material  will  become  flatter,  while  it  will  be 
steeper  in  the  less  permeable  material  hence  inducing  higher 
seepage  forces . 

No  explanation  has  been  given  of  the  fact  that  during 
the  first  filling  the  movements  started  at  reservoir  eleva¬ 
tion  635  m.  and  were  very  high  at  650  m.  elevation  but  that 
during  the  second  filling  the  movements  started  at  reservoir 
elevation  650  m.  only.  If  the  hypothesis  of  constant  dis¬ 
charge  is  accepted  it  can  be  argued  that  since  the  filling 
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of  the  reservoir  results  in  the  raising  of  the  water  table, 
the  movements  should  have  started  again  at  reservoir  eleva¬ 
tion  635  m.  during  the  second  stage  of  filling.  This  is 
explained  by  the  simple  fact  that  the  first  reservoir  fill¬ 
ing  up  to  635  m.  and  then  to  the  650  m.  elevation  was  done 
during  the  wet  period  of  the  year  (October,  1960  -  November, 
1960).  The  precipitation  is  generally  high  during  these 
months  in  that  area  and  as  the  recharge  of  the  ground  water 
basin  will  be  more,  it  is  obvious  that  the  discharge  passing 
out  will  be  more  as  compared  to  the  other  periods  of  the 
year.  This  means  that  during  this  period  of  the  first  fill¬ 
ing,  the  mountain  water  table  was  slightly  higher  and  steeper, 
passing  slightly  higher  discharges  than  in  the  normal  periods. 
During  the  second  stage  of  filling  the  reservoir  level  of 
635  m.  was  reached  during  the  colder  month  of  January,  1961, 
which  is  a  comparatively  dry  period  and  thus  it  is  expected 
that  during  this  period  the  mountain  water  table  must  be 
slightly  lower  than  the  previous  case  as  the  discharge  pass¬ 
ing  out  is  less.  Therefore,  it  is  evident  that  at  the  re¬ 
servoir  elevation  of  635  m.  the  same  conditions  did  not 
occur  again,  the  water  table  is  still  lower  than  the  pre¬ 
vious  case  and  so  the  factor  of  safety  is  also  slightly 
more.  The  factor  of  safety  became  close  to  unity  again 
when  the  reservoir  elevation  was  at  650  m.  and  thus  the 
movements  started  again.  The  same  phenomenon  occurred  again 
in  1963  during  the  second  lowering  and  third  filling. 

Moreover,  as  the  permeability  of  the  material  which  was 
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submerged  earlier  has  increased  the  increase  in  head  in 
the  material  of  two  different  permeabilities  material 
will  not  be  the  same  as  in  the  previous  case.  The  slope 
of  the  mountain  water  table  will  become  flatter  if  the  per¬ 
meability  has  increased,  which  means  less  seepage  forces 
than  the  previous  case. 

Muller  (1968)  has  also  mentioned  that  although  heavy 
preci pi tati on  occurred  in  the  months  of  August  and  September 
1963  prior  to  failure,  no  correlation  with  the  piezometric 
observations  could  be  established.  This  shows  that  some 
distinct  changes  in  the  permeability  of  the  rock  material 
along  the  inclined  and  horizontal  beddings  might  have  taken 
place. 

Another  question  still  remains  unresolved.  The  pie¬ 
zometric  observations  during  the  later  part  of  filling  do 
not  show  any  correlation  between  a  raising  of  the  mountain 
water  table  with  the  raising  of  the  lake  level. 

This  may  be  explained  by  the  fact  that  during  the 
first  filling,  no  piezometric  observations  were  made,  as 
such  the  initial  slope  of  the  mountain  water  table  was  not 
known.  As  the  submergence  proceeded,  the  permeability  of 
the  submerged  rock  material  increased  considerably,  while 
the  inclined  layers  of  Malm  and  lower  cretaceous  still  had 
low  permeab i 1 i ty .  It  therefore  seems  necessary  that  the 
mountain  water  table  should  have  a  steeper  gradient  in  these 
layers  than  the  ones  which  were  submerged.  No  data  of  the 
permeability  as  well  as  the  mountain  water  table  in  the  in- 
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dined  layers  of  Malm  and  lower  cretaceous  is  available. 

However,  it  appears  from  the  locations  of  these  pie¬ 
zometers  that  they  were  observing  the  water  levels  in  the 
more  permeable  material  where  the  bedding  was  almost  hori¬ 
zontal.  Hence  it  is  evident  that  these  are  not  representa¬ 
tive  observations  of  the  mountain  water  table  in  the  less 
permeable  inclined  layers  of  Malm  and  lower  cretaceous 
formation. 

A  hypothetical  free  surface  line  in  these  layers  is 
shown  on  Fig.  5.1,  The  x-section  is  taken  from  Muller's 
report  (1968)  and  gives  the  possible  explanation  for  not 
observing  higher  mountain  water  table  in  this  region.  As 
shown  the  gradient  of  the  water  table  in  the  inclinded 
layers  is  steep,  resulting  in  high  pore  pressures  and  seep¬ 
age  forces,  while  the  water  table  in  the  horizontal  beds 
where  the  piezometers  were  installed  is  quite  flat. 

It  is  also  interesting  to  note  that  depending  upon 
the  failure  profile  chosen,  several  investigators  have  found 
the  values  of  angle  of  shearing  resistance  at  failure  rang¬ 
ing  from  17.5°  to  24°  (Muller  (1968)).  Nonveiller  (1967) 
has  also  found  the  values  of  angle  of  shearing  resistance 
at  failure  from  27°  to  28.5°,  but  the  assumptions  of  slip 
surface  position  and  form  differ  very  much  from  nature 
(Muller  (1968)).  Kenney  (1967)  obtained  a  value  of  angle 
of  shearing  resistance  less  than  22°  at  failure.  These 
values  are  much  smaller  than  what  could  be  expected  for  lime¬ 
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FIG.  5.1  DIAGRAM  SHOWING  GEOLOGICAL  SECTION  OF  THE  VAJONT  ROCK  SLIDE  WITH  HYPOTHETICAL  WATER  TABLE 

(GEOLOGICAL  SECTION  TAKEN  FROM  FIG.  2.  MULLER(1968) ) 


93 


According  to  Skempton  (1966)  a  friction  angle  for 
limestone  of  about  30°  should  be  expected.  He  is  also  of 
the  opinion  that  in  the  case  of  Vajont  a  reduction  of  fric¬ 
tion,  as  a  result  of  slip  between  adjacent  beds  should  not 
be  considered,  as  none  of  the  observations  showed  any 
smoothing  of  the  roughness  pattern  on  the  bedding  plane. 

It  is  now  evident  that  these  low  values  of  friction 
which  are  not  representati ve  values  for  the  limestone  are 
obtained  because  in  all  the  stability  analyses,  the  water 
table  is  assumed  more  or  less  horizontal  and  equal  to  the 
reservoi r  1 e vel  . 

If  high  seepage  forces  as  a  result  of  a  steeper  water 

table  in  the  inclined  layers  of  Malm  and  lower  cretaceous 

formations  are  taken  into  account,  it  seems  necessary  to 
obtain  a  value  of  angle  of  shearing  resistance  close  to  the 
representative  value  of  limestone. 

To  check  this,  an  approximate  stability  analysis  based 
on  Janbu's  ejt  aj_.  (  1  956  )  method  is  done.  The  valley  cross- 
section  and  the  slip  surface  chosen  correspond  to  slip  sur¬ 
face  II,  Fig.  2  in  Kenney's  paper.  The  pore  pressure  in¬ 
serted  in  the  stability  analysis  are  based  on  the  hypo¬ 
thetical  free  surface  line  shown  in  Fig.  5.1.  The  stability 

calculations  indicate  that  when  the  rock  slope  failed  the 

average  coefficient  of  sliding  resistance  was  0.554..  This 
value  correspond  to  an  angle  of  shearing  resistance  of  29°, 
which  is  obviously  a  more  representative  value  for  lime¬ 


stone. 


fl 
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To  summarize  it  is  concluded  that  this  possible 
explanation  of  the  Vajont  rock  slide,  based  on  the  constant 
discharge  concept,  provides  a  qualitative  answer  to  most 
of  the  unresolved  questions  pertaining  to  the  complexities 
of  this  major  catastrophe.  At  the  same  time  it  is  rea¬ 
lized  that  in  the  absence  of  any  data  on  the  permeabilities 
of  different  layers  and  the  mountain  water  table  in  the  in¬ 
clined  layers,  this  explanation  is  purely  hypothetical. 
However,  it  may  contribute  to  an  understanding  of  the  im¬ 
portance  of  the  detailed  knowledge  of  ground  water  hydrology 
and  in-situ  permeability  in  natural  slope  stability  problems. 


■ 


LIST  OF  REFERENCES 


Allen,  D .  ,  1954.  "Relaxation  Method".  McGraw-Hill  Book 
Co .  ,  Inc  .  ,  New  York. 

Aravin,  V.I.  and  Numerov,  S.N.,  1965.  "Theory  of  fluid 

Flow  in  Unde formab le  Porous  Media".  Israel  Program 
for  Scientific  Translations,  Jerusalem. 

Bear,  J.?  Zaslavsky,  D.  and  Irmay,  S.,  1968.  "Physical 

Principles  of  Water  Percolation  and  Seepage".  UNESCO, 
Pari  s  . 

Bishop,  A.W.,  1952.  "The  Stability  of  Earth  Dams".  Ph.D. 
Thesis,  University  of  London. 

Bishop,  A.W.,  1955.  "The  Use  of  the  Slip  Circle  in  the 
Stability  Analysis  of  Slopes".  Geo  technique,  5, 

7-17. 


Bishop,  A.W.  and  Bjerrum,  L.,  1960.  "The  Relevance  of  the 
Triaxial  Test  to  the  Solution  of  Stability  Problems". 
Proc.  Research  Conf.  Shear  Strength  of  Cohesive  Soils, 

(Boul  der-,  Col  arado)  ,  437-501  . 

Casagrande,  L.,  1932.  "Naeherungsmethoden  zur  Bestimmung 
von  Art  und  Menge  der  Sickerung  durch  geschuettete 
Daemme".  Thesis,  Technische  Hochschule,  Vienna. 
Translated  by  U.S.  Corps  of  Engineers,  Waterways, 

Exp.  St a.,  Vicksburg,  Miss. 


Charmonman,  S.,  1967.  "Coastal  Parallel  Canals  with  Inter¬ 
mediate  Drains".  J.  Hydraul .  Div.,  Proc.  ASCE,  93(HY-1), 
13-27. 


Cheung,  Y.K.,  Watson,  M.  and  Skjolingstad,  L.,  1970. 

"Steady  State  Seepage  and  Transient  Ground  Water  Flow 
in  Slopes".  23rd  Annual  Canadian  Geotechnical  Con- 
fe ren ce ,  Banff,  Alberta. 

DeLory,  F.A.,  1957.  "Long  Term  Stability  of  Slopes  in 

Over-Consolidated  Clays".  Ph.D.  Thesis,  University 
of  London . 

Dupuit,  J.,  1863.  "Etudes  Theoriques  et  Pratiques  Sur  le 
mouvement  des  eaux  dans  les  Canaux  decouverts  et  a 
travers  les  terrains  Permeable".  Paris. 


96 


Finnemore,  E.J.  and  Perry,  B.,  1968.  "Seepage  through  an 
Earth  Dam  Computed  by  the  Relaxation  Technique". 

Water  Resources  Research,  4(5),  1059,  Washington,  D.C. 

Finn,  W.D.  Liam,  1967.  "Finite  Element  Analysis  of  Seepage 
through  Dams".  J.  Soil  Mechs.  and  Found.  Div.  Proc. 
ASCE  ,  93,  SM6,  4 T.  —  ' 


Forsythe,  G.E.  and  Wasow,  R.W.,  1960.  "Finite  Difference 
Methods  for  Partial  Differential  Equations".  John 
Wiley  and  Sons  ,  Inc .  ,  New  York. 

Harr,  M.E.,  1962.  "Ground  Water  and  Seepage".  McGraw-Hill 
Book  Co . ,  Inc.,  New  York. 

Henkel,  D.J.  and  Skempton,  A.W.,  1955.  "A  Landslide  at 
Jackfield,  Shrophire,  in  Heavily  Over-Consolidated 
Clay".  Geotechni  que  ,  5,  131. 

Jaeger,  C.  ,  1  969.  "The  Stability  of  Partly  Immerged  Fissured 
Rock  Masses,  and  the  Vajont  Rock  Slide".  Civ i 1 
Engineering  and  Public  Work  Review,  1204-1207. 

Janbu,  N.,  Bjerrum,  L.  and  Kjaernsli,  B.,  1956.  "Veiledning 
Ved  losning  av  f un damen te r i n gs o pp ga ve r " .  Norwegi an 
Geotechnical  Institute  Publication  No.  16,  Oslo. 

Jeppson,  R.W.,  1968.  "Seepage  through  Dams  in  the  Complex 
Potential  Plane".  J.  Irrig.  Drainage  Div.,  Proc. 

ASCE,  V  ,  94 ( I RI ) ,  23-39. 

Jeppson,  R.W.,  1968.  "Seepage  from  Ditches  -  Solution  by 
Finite  Differences".  J.  Hydraul  .  Div.,  Proc.  ASCE, 

V  .  94 ( HY -1  ) ,  259-283. 

Kenney,  T.C.,  1967.  "Stability  of  the  Vajont  Valley  Slope". 
Rock  Mechanics  and  Engineering  Geology,  Vol .  VII, 

10-16. 

Linsley,  Jr.,  R.K.,  Kohler,  M.A.  and  Paulhus,  Joseph  L.H., 
1958.  "Hydrology  for  Engineers".  McGraw-Hill  Book 
Co .  Inc.,  New  York. 

McNown,  J.S.,  Hsu,  E.Y.  and  Yih,  C.S.,  1955.  "Application 
of  the  Relaxation  Technique  in  Fluid  Mechanics  with 
Discussion  by  Others".  T rans  c .  ASCE,  Vol .  1  20,  650-686. 

Morgenstern,  N.R.  and  Price,  V.E.,  1965.  "The  Analysis  of 

the  Stability  of  General  Slip  Surfaces".  Geotechni que , 
15,  79-93. 


97 


M.I.T.,  1969,  Department  of  Civil  Engineering,  Soil  Mechani cs 
Publication  No.  235. 

Muller,  L.,  1964.  "The  Rock  Slide  in  the  Vajont  Valley". 

Rock  Mechanics  and  Engineering  Geology,  Vol  .  II, 
148-212. 

Muller,  L.,  1968.  "New  Consideration  on  the  Vajont  Slide". 
Rock  Mechanics  and  Engineering  Geology,  Vol.  VI, 

1-91  . 

Neuman,  S.P.  and  Witherspoon,  P.A.,  1970.  "Finite  Element 

Method  of  Analyzing  Steady  Seepage  with  Free  Surface". 
Water  Resources  Research,  Vol.  6,  No.  3,  889-897. 

Nonveiller,  E.,  1967.  "Zur  Frage  der  Fel s tru ts chung  in 
Vaiont-Tal".  Fel sm .  u .  Ing  .  Geo!  .  ,  Vol.  V / 1  . 

Pavlovsky,  N.N.,  1931.  "Seepage  through  Earth  Dams". 

Instit.  Gi drotekhni ki  i  Mehoratsii,  Leningrad. 
Translated  by  U.  S'.  Corps  of  Engineers . 

Pol ubari no va-Kochi na  ,  P.Ya.,  1  962.  "Theory  of  Ground  Water 
Movement".  Translated  by  J.M.  DeWiest,  Princeton, 

N.J. 

•  • 

Schaffernak,  F.,  1917.  "Uber  die  s tands i cherhei t  durch- 
laessiger  geschuette ter  Damme".  Allgem.  Bauzeitung. 

Shaw,  F.S.  and  Southwell,  R.V.,  1941.  "Relaxation  Methods 

Applied  to  Engineering  Problems".  Proc.  Royal  Society 
London  ,  London. 

Skempton,  A.W.,  1964.  "Long  Term  Stability  of  Clay  Slopes". 
Geotechni gue ,  1 4 ,  77-101. 

Skempton,  A.W.,  1966.  "Bedding  Plane  Slip,  Residual  Strength 
and  the  Vajont  Landslide".  Geo te cnn i que  ,  16,  82-84. 


Taylor,  R.L.  and  Brown,  C.B.,  1967.  "Darcy  Flow  Solutions 
with  Free  Surface".  J.  Hydraul .  Div.,  Proc.  ASCE, 
Vol .  93 (HY-2 )  ,  25-33. 

Thom,  A.  and  Apelt,  C.J.,  1961.  "Field  Computations  in 

Engineering  and  Physics".  D.  Van  Nostrand  Co .  Ltd. , 
London. 

Todd,  D.K.,  1959.  "Ground  Water  Hydrology".  John  Wiley 
and  Sons  Inc.,  New  York. 


98 


Verruijt,A.,  1970.  "Theory  of  Ground  Water  Flow". 

Macmillan,  London. 

Zienkiewicz,  0.  and  Cheung,  Y.K.,  1965.  "Finite  Element 
in  the  Solution  of  Field  Problems".  The  Engineer, 
London,  280. 

Zienkiewicz,  0.,  Mayer,  P.  and  Cheung,  Y.K.,  1966. 

"Solution  of  Anisotropic  Seepage  by  Finite  Elements". 
J.  Eng.  Mech.  Div.,  ASCE,  92,  EMI,  111-120. 


' 


APPENDIX  A 


NOMENCLATURE  FOR  COMPUTER  PROGRAM 
-  LISTING  OF  THE  COMPUTER  PROGRAM 


A- 1 


H  (  ) 
Hs(  ) 


X(  ) 

V(  ) 
RES  (  ) 
CH K(  ) 


KGP  (  ) 


NPP (  ) 

N 

Q 

DELTA 

rMAX 

JMAX 

ME 

ARESM 


Nomenclature  for  Computer  Program 

Grid  point  potential. 

Potential  at  the  free  surface;  2nd  and 
3rd  column  of  the  same  matrix  stores  inter¬ 
polated  values  of  potentials  in  the  vertical 
direction  corresponding  to  a  particular 
point  on  the  free  surface. 

Abscissa  of  the  grid  point. 

Ordinate  of  the  grid  point. 

Residue  at  each  grid  point. 

Matrix  which  ensures  convergence  of  the 
free  surface  potentials  within  the  allow¬ 
able  limits  of  error. 

Array  which  stores  the  specific  free  surface 
points  which  are  coincident  with  the  grid 
points. 

Array  which  stores  the  actual  grid  points 
which  may  correspond  to  the  free  surface  points, 
Number  of  free  surface  points. 

Dimensionless  infiltration  factor. 

Mesh  size. 

Maximum  number  of  rows  in  the  H(  )  matrix. 
Maximum  number  of  column  in  the  H(  )  matrix. 
Number  of  linear  equations  to  be  solved  by 
subroutine  SOLVE. 

Maximum  allowable  residue. 


I 


MI  TR 

EPS  LN 
ALERO 
RELAX 

PINT 

SOLVE 

ARRA 

SIMQ 

LCNT 

KGPCH 

LUT 

KAN  (  ) 

YY 

KQQ 

KA 


A-2 


Maximum  number  of  iterations  required  for 
convergence . 

Over-rel axa ti on  factor. 

Allowable  error  in  the  free  surface  potentials. 
Name  of  the  subroutine  used  for  relaxation 
by  SOR. 

Name  of  the  subroutine  used  for  interpola¬ 
tion  of  grid  potentials. 

Name  of  the  subroutine  used  for  assigning 
values  of  coefficients  and  constants  to 
be  used  in  subroutines  ARRA  and  SIMQ. 

Name  of  the  subroutine  which  converts  a 
matrix  into  an  array  (SSP). 

Name  of  the  subroutine  used  for  solving 
linear  equations  simultaneously  (SSP). 
Instantaneous  count  for  the  number  of  itera¬ 
tions  in  subroutine  RELAX  before  achieving 
convergence . 

Index  used  to  trace  the  points  in  the  array 
KGP . 

Index  for  storing  points  for  array  KGP. 

Array  used  for  deciphering  i nterpol ati on 
in  X  and  Y  directions. 

Variable  used  to  reduce  the  ordinate  values 
in  steps  of  delta. 

Index  for  array  KGP. 

Index  for  array  KAN  .. 


V 


A-3 


NL1 
ML  1 
A(  ) 

C(  ) 


ALPHA 


DELT-j  ,  DELT^ 


delt3>  delt4 


JAR-j  ,  J AR^ 


KAJ 

C  ( 5 ) 


Index  used  for  a  particular  row  in  matrix 

H (  ). 

Index  used  for  a  particular  column  in 
ma tri x  Hs (  ) . 

Matrix  for  coefficients  used  in  the  linear 
equations. 

Array  for  storing  constants  of  simultaneous 
linear  equations. 

Angle  subtended  by  the  segment  of  the  free 
surface  line  with  the  horizontal. 

Horizontal  distances  of  grid  points  which 
are  interpolated  or  otherwise  used  in  the 
polynomial. 

Vertical  distances  of  the  grid  points  which 
are  interpolated  or  otherwise  used  in  the 
polynomi al . 

Indices  used  for  those  surface  potentials 
that  are  required  to  be  used  as  constant 
in  the  solution  of  polynomial. 

Index  used  for  specifying  column  number  in 
the  ma  tri x  H (  ) . 

Constant  for  infiltration  factor. 


o  n 
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MA  I  N 


MA  I  N 


C  INFILTRATION 

COMMON  H(30>60) >HS(60t3) »  X ( 60 )  *Y(60)  *  RES  C  500 )  >CHK(60 
*•160)  >  KGP (25)  , 

1NPP ( 60 ) 

COMMON  N  » Q  >  DELTA 
C  INPUT  DATA 

READ  (5*10) I  MAX  >  JMAX * N  * NE  »  DELTA  > Q 

10  FORMAT ( 4  I  4  >  2F10 . 3  ) 

WRITE (6  >137) I  MAX > JMAX » N >NE > DELTA  ,0 

137  FORMAT ( ,1*>////>4I5>F10.4>F10.4) 

READ ( 5  > 11 ) ARESM  >  M I TR  *  EPSLN  * ALERO 

11  FORMAT  ( F10  a  3  > I  4  >  2F10 . 3 ) 

WR'ITE(6  >133  )ARESM»MITR»EPSLN*ALERO 

138  FORMAT ( '0'>F10.4>I5>2F10.4) 

READ  (5  >12 )  (X ( I  )  >  1  =  1 >N ) 

12  FORMAT  ( 8 F 1 0 • 4 ) 

WRITE (6  >145 )  (X ( I  )> I  =  1 >N) 

145  FORMAT (*0’  >4E13a3  ) 

READ  (5  >12 )  ( Y ( I  )  > 1  =  1 >N) 

W  R  I  T  E  (  6  >  1 4  5  )  (  Y  (  I  )  >  I  =  1  >  N  > 

INITIALIZATION  OF  FREE  SURFACE  AND  INTERNAL  GRID 
POTENTIALS 

J  =  1 

DO  20  I  =  1 > I  MAX 
H( I » J )  =  1400. 

20  CONTINUE 

DO  31  I  =  1»N 
HS (  I  *  1 )  =  Y(I) 

31  CONTINUE 

DO  25  I  = 1 > I  MAX 
DO  25  J  =  2  >  JMAX 
H ( I >  J  )  =  500. 

25  CONTINUE 
NM  =  N-28 

DO  26  I  =NM*N 
HS ( I > 1  )  =  1000. 

26  CONTINUE 
H ( I  MAX  •  JMAX )  =  1000. 

WRITE  ( 6 > 30 ) ( HS ( I > 1 ) > I =1 » N ) 

30  FORMAT  ( 'O’ >4E13.3) 

JCHK  =  0 

777  JCHK  =  JCHK  +  1 

IF (JCHK-160 )991  >992  >992 
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MAIN  ...  (CONT’D) 


991  DO  75  5  JMP  =  1*N 

CHK ( JMP  » JCHK ) =  HS ( JMP > 1 ) 

755  CONTINUE 

IF(JCHK-l)2f754»2 
2  ERMAX  =  0.0 

DO  733  JERR  =  12*44 

ERR1  =  ABS(CHK(JERR*JCHK)-CHK(JERR*JCHK-1) ) 

IF ( ERR1  -  ERMAX ) 733»711 *711 
711  ERMAX  =  ERR  1 
733  CONTINUE 

IF (ERMAX-ALERO) 7  00  >700  >7  54 
700  I CHK  =  0 

GO  TO  789 

754  CALL  RELAX ( I  MAX  *  JMAX  »  EPSLN  > ARESM  *  Ml TR ) 

CALL  PINT ( I  MAX  »  UMAX ) 

CALL  SOLVE(NE*IMAX*jMAX) 

GO  TO  777 
789  WRITE  (6*100) 

100  FORMAT* ' 1 ' *10X* ’ POTENTIALS  AT  GRID  POINTS’*///) 

DO  103  1  =  1*1  MAX 

WRITE  (6*101)  (H( I »J)  »  J= 1 »  JMAX ) 

101  FORMAT  (8 (F10.4.5X) ) 

103  CONTINUE 

WRITE  (6*105) 

105  FORMAT  (’ 0 ’* 20X »’ POTENT  I ALS  AT  THE  FREE  SURFACE’*////) 
WRITE  (6*106) (HS( I *1) >1=1 >N) 

106  FORMAT  (8(F10.4*5X) ) 

WRITE  (6*107) 

107  FORMAT(  *  1  *  »20X> 'COORDINATES  OF  THE  POINTS  ON  FREE 
*  SURFACE’ *///) 

DO  108  1  =  1  *  N 

WRITE  (6*109)  V ( I ) >X( I ) 

1 09  FORMAT  ( ’ 0 ’ * 10X » F10 .4 , 5X * F10.4 ) 

108  CONTINUE 

992  WR I TE ( 6  * 993 ) 

993  FORMAT ( ’0’ *10X* ’ FAILED  TO  CONVERGE  IN160  ITERATIONS') 
STOP 
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SUBROUTINE  RELAX 


SUBROUTINE  RELAX  (  I  MAX  *  JMAX  »  EPSLN  * ARESM  *  Y  I  TR  ) 

COMMON  H<30  *60)  »HS( 60  *3)  *X(  60) >Y( 60) *  RES (500) »CHK(60 
**160)  *  KGP ( 2  5 )  , 

1 v  P  P  (  6  0  ) 

COMMON  N  *  Q  * DELT A 

THIS  SUBROUTINE  CALCULATES  POTENTIALS  BV  USING  SOR 

LUT  =  1 

DO  250  IP  =  1*25 
KGP  ( IP )  =  0 
250  CONTINUE 
NPP(l)  =  0 
DO  13  LUP  =  2 » N 
XK  =  X(LUP) 

9  ALAX  =  XK/DELTA 

I F ( AL AX  -  1. ) 13*14*16 
16  XK  *  XK  -  DELTA 
GO  TO  9 
14  YJ  =  Y(LUP) 

10  ALAY  =  YJ/DELTA 

I F ( ALAY  -  1 «  )  13  »  3  9  *  40 

39  KGP(LUT)  =  LUP 
LUT  =  LUT  +  1 
GO  TO  13 

40  YJ  =  YJ  -  DELTA 
GO  TO  10 

13  CONTINUE 
M  =  2 
YY  =  Y ( 1 ) 

LUT  =  1 
DO  50  L  =  2 » N 

IF  { ( YY-Y ( L ) ) -DELTA )  51>52»52 

51  NPP(M)  =  L 

M  -y+1 

GO  TO  50 

52  YY  =  YY-DELTA 

IF  (L-KGP(LUT)  )  50  *  53  »  5  0 

53  NPP(M)  *  L 
LUT  =  LUT  +  1 
M :  =  M  +  1 

50  CONTINUE 
KGPCH  =1 
LCNT  =  0 
222  HT  =  0,0 

LCNT  =  LCNT+1 
YY  =  Y ( 1 ) 

J  =  1 

27  J  =  J+ 1 
I  =  I  MAX 
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SUBROUTINE  RELAX  ...(CONT'D) 

93  IF ( J-( JMAX-1) )54»261*260 
261  NM  =  N-l 

AN  =  (Y(NM)-Y(N) ) /DELTA 

HO  =  ( HS ( NM %  1 ) / ( 1 » +AN##2 )  )  +  ( A  N *  *  2 / ( 2 • * ( 1 • +AN*#2 )  )  ) 

** ( H ( I  MAX  >  JMAX-2 ) 

1+H (  I  MAX  » JMAX )  ) 

DH  =  ABS ( HQ-H ( IMAX  * JMAX-1 )  ) 

IF  (HT-DH) 88  *  37 >87 
88  HT  =  DH 

87  H ( IMAX > JMAX-1 )  =  H ( IMAX ♦ jMAX-1 ) +EPSLN* ( HO-H ( I  MAX > JMAX 
*-l  )  ) 

GO  TO  27 
54  XX  =  DELTA 

C  RELAX  ALONG  BASE  LINE 

HO  =  ( H ( I  MAX *  J  +  l )+2.*H( IMAX-1 »j)+H( IMAX* J-l)  )*.25 
DH  =  ABS ( HO-H ( IMAX*j) ) 

IF (HT-DH) 231  *232  *232 
231  HT=DH 

2  32  H(  I  MAX  *  J ) =H ( IMAX » J ) +EPSLN# ( HO-H ( IMAX*J)  ) 

25  I  =  I  -  1 
47  XX  =  XX+DELTA 
KJ  =  NPP(J) 

IF  (XX-Y(KJ)  ) 20  »  2 1 » 21 

C  RELAX  INTERIOR  POINTS 

2  0  HO= ( H ( I  * J  +  l )+H( 1-1 » J)+H( I » J-l )+H ( 1  +  1  * J ) )*0.25 
DH=ABS ( HO-H ( I »j) ) 

IF(HT-DH)215*214*214 
215  HT  =  DH 

214  H  (  I  *  J  )  =H  (  I  *  J)+EPSL»N*(HO-H  (  I  »  J  )  ) 

GO  TO  25 
21  KJ  =  NPP(J) 

IF  (  ( X ( KJ+1 ) -X ( KJ )  )“DELTA)67  *30*48 
67  IF  (KJ-KGP(KGPCH) )31*29»31 
31  EM DA  =  (X(KJ+1)-X(KJ) ) /DELTA 
AN  =  (DELTA-(XX-Y(KJ) ) )/DELTA 

HO= ( ( 2 • *HS ( KJ+1 ♦ 1 ) ) / ( EMDA* ( 1 • +EMDA ) )  +  ( 2 . *HS ( K J *  1 ) ) / ( AN 
*  *  ( 1  •  +  A  N  )  )  +  (  2 

l.#H(  I  *J“1  )  )  /(  l.+EMDAK  (  (2*#H(  1+1  »j>  )  /  (  l.+AN)  )  )  *  (  EMDA 

#*AN)/(2.#(EMDA 

1+AN ) ) 

DH  =  ABS ( HO-H ( I » J ) ) 

IF  ( HT-DH ) 18  » 19  *  IS 

18  HT  =  DH 

19  H ( I  *  j )  =  H( I *J)+EPSLN*(HO-H( I  ,j)  ) 

GO  TO  27 

29  EMDA  =  (X(KJ+1)-X(KJ) ) /DELTA 
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SUBROUTINE  RELAX  ...(CONT'D) 

HO- ( ( 2 • *HS ( K  J  +  l * 1 )  ) / ( EMDA* ( 1 • +  EMDA ) >  +  ( 2 • *H ( I  * J-l )  > / ( 1 • 
*+E^DA)+HS(KJ 

1»1)+H(I+1»J)  )*(  EM.DA  /  (  2  •  #EMDA  +  2  •  )  ) 

DH  =  ASS ( HO-H (I  »  J  )  ) 

IF  ( HT-DH ) 4 1 » 42  » 42 

41  HT=DH 

42  H  (  I  »  J  )  =  H( I » J ) +EP5LN* ( HO-H (  I  » J  )  ) 

KGPCH  =  KGPCH+1 

GO  TO  27 

30  AN  =  (DELTA-(XX-Y(KJ) )) /DELTA 

HO=( (  (2.#HS(KJ*1)  ) /(AN*( l.+AN)  ) )  +  (  (  2  •  *  H  ( I  + 1  *  j ) )/(le 
*+ AN )  ) +H (  I » J+l ) + 

1H (  I »  J-l )  ) * ( AN/ ( 2 •* ( AN+1 • )  )  ) 

DH  =  ABS ( HO-H ( I #J) ) 

IF  ( HT-DH ) 99  » 22  *  22 
99  HT  =  DH 

22  H ( I » J  )  =  H ( I > J ) +EPSLN* ( HO-H ( I  * j )  ) 

GO  TO  27 

260  IF  (LCNT-MITR)270»270»273 

270  RES ( LCNT ) =HT 

IF  (HT-ARESM)271>271,222 

271  WRITE (6  *37 )  RES ( LCNT )> LCNT 

37  FORMAT  ( ' 0 ' » 20X » 1  MAX  RESIDUE^'  » E 12  *  5  »  2X  » '  ITERATION 
*  N0= '  » 14 ) 

GO  TO  273 
48  WRITE  ( 6  ? 60 ) 

60  FORMAT (' 1 '* 10X »' UNEXPECTED  HAPPENEDIN  RELAX') 

273  RETURN 
END 
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SUBROUTINE  PINT 
SUBROUTINE  P I  NT ( I  MAX >  JMAX ) 

C  THIS  SUBROUTINE  INTERPOLATES  THE  VALUES  OF  POTENTIALS 

DIMENSION  KAN ( 20 ) 

COMMON  H(30>60)  *  HS (60*3)  >X(60) »Y(60)  »RES(500)  >CHK(60 
**160)  »  KGP ( 2  5 )  , 

1NPP ( 60 ) 

COMMON  N»Q» DELTA 
C  INITIALIZATION 

1  =  2 

NM3  =  N-29 
KQG  =  1 

YY  =  Y ( 1 ) 

J  =  1 
JJ  =  1 
KA  =  J  J 
JL  =  1 

HS ( 1 »  2  )  =  H ( 1 ?  1  ) 

HS (1*3)  =  H ( 1  *  1  ) 

DO  30  LP  =  1*20 
KAN(LP)  =  0 

30  CONTINUE 
20  J=J+1 

JL  =  JL  +  1 
K  J  =  J 

I F ( j~ ( N~28  )  )  15  *  100*  100 
15  I F (  ( YY-Y ( J )  )  -  DELTA ) 1  ♦  1 »  3 
3  YY=YY-DELTA 
JL  =  JL  -  1 

I F ( J  -  KGP(KQO)  ) 3 1  *  32  *  3 1 
32  HS ( j  »  2 )  =  H ( I  +  1  *  JL ) 

HS ( J  »  3 )  =  H ( 1+2  » JL) 

KOQ  =  KQQ  +  1 
GO  TO  20 

31  WRITE (6 *51 )X( J)  »X(j-2)  *  DELTA 

51  FORMAT (  *  0 •  *4X* '  IN  P  I  NT  '  » 3 F 1 3 . 5 ) 

I F ( (X(j)-X(J-2) ) -DELTA) 15 *5 *50 
1  IF( (X(j)-X(J-l) )-DELTA)4*5»50 

C  INTERPOLATION  IN  Y-DIRECTIONS 

5  IF(KAN(KA)-J)25 *17*25 
17  KJ  =  KJ  +  1 

25  DO  6  11=1*2 

NL  =  KJ  +  II 

ML  =  II  +  1 

HS(J>  ML  )  =  HS(J»1)  “  (  ( HS ( J  » 1 ) -H ( I  *  J  L )  ) * ( Y ( J ) -Y (  NL 
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SUBROUTINE  PINT  ..*(CONT»D) 

#  ))/( DELTA- 
1 ( YY-Y ( J )  )  )  ) 

36  I  F  <  (YY-Y(KJ  +  2) )-DELTA)6*7*23 
7  KAN(JJ)  =  J 
JJ= JJ  +  1 
!<  A  =  J  J 

KAN(JJ)  =  J+1 
J  J  =  J  J  + 1 
K J  =  J+ 1 

23  ml  =  ML  +  1 
IK  =  2 

27  I F ( J-NPP ( IK) ) 18  *  1 9  *  1 8 
13  IK  =  IK  +  1 

GO  TO  27 

19  IF  ( (YY-DELTA)-{ Y(Kj+2)+DELTA)  ) 28  O  52  *  28 

28  XY  =  (  ( YY-DELTA ) -Y ( K J+2 )  ) 

GO  TO  29 

29  HS ( J  *  ML )  =  H( I  *  IK )-( (H( I  * IK)-H( 1+1  ♦  IK)  >*XY) /DELTA 
GO  TO  20 

6  CONTINUE 
GO  TO  20 

152  HS( J»YL)  =  H ( I  + 1  » I  K  ) 

GO  TO  20 

INTERPOLATION  IN  X-DIRECTION 

4  DO  10  12=1*2 
NL1  =12+1 
ML1  =  12  +  1 

HS ( J  *  ML  1 )  =  H (  NL1*JL-1)  -  (  ( H ( N  L 1  *JL-1)-H(  NL1*JL)) 
#*( DELTA- ( X (K 

1J+1)  -  X ( K J ) ) ) / ( X ( K J+1 5  -  X ( K J-l ) ) 5 
10  CONTINUE 
I  =  1+1 
GO  TO  20 
50  WRITE  (6,60) 

60  FORMAT! *1* »10X» ‘DIFFERENCE  IN  X  EXCEEDS  DELTA  IN  RINT 

*  SUBROUTINE  *  ) 

100  RETURN 

END 
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SUBROUTINE  SOLVE 

SUBROUTINE  SOLVE ( NE * I  MAX *  JMAX ) 

DIMENSION  A ( 5  *  5 )  *  S  A  (  2  5  )  *C(5) 

COMMON  H  (  3  0  *  6  0  )  »  H  S  (  6  0  »  3  )  »  X ( 6  0  )  *  Y  {  6  Q  )  »RES(500)  ,CHK(60 
**160)  *  KGP ( 2  B ) * 

1NPP(60) 

COMMON  N*Q* DELTA 


THIS  SUBROUTINE  SOLVES  THE  LINEAR  EQUATIONS 


DO  10  1  =  1 »  N E 
DO  10  J  =  1  *  N E 
A  (  I  ♦  J  )  =  0 . 0 
10  CONTINUE 
KQA  =  1 
1  =  1 


KJ  =  2 
JL  =  2 
YY  =  Y ( 1 ) 

NM1  =  N-28 
M  =50 

A  (  1  *  1  )  =  1 . 0 
A  (  1 » 4  )  =  0  •  0 
A  ( 1  *  5 ) =0 • 0 
A  (  2  *  1  )  =  1  •  0 
A  (  2  *  4  )  =  0  o  0 
A  (  2  *  5  )  =  0  •  0 
A  (  3  » 1  )  =  1 8  0 
A(3*2 )=0.0 
A (3*3) =0.0 
A  (  4  » 1  )  =  1 . 0 
A  (  4  *  2  )  =  C  e  0 
A  (  4  *  3  )  =  0  8  0 
A  (  5  *  1  )  =  0  8  0 
A  (  5  *  3  )  =  0 . 0 
A(5*5)=080 
DO  10  0  J  =  2  *  NM 1 

ALPHA=  ATAN( (Y( j-1)-Y( j) ) / ( X ( J)-X( j-1) ) ) 
A ( 5  *  2 ) =S I N ( ALPHA  ) 

A ( 5  *4 ) =COS ( ALPHA ) 


IF( J-2)71, 71*77 
71  DELT1  =  DELTA 
DELT2  =  DELTA 
DELT3  =  DELTA-  ( YY  -  Y ( J ) ) 
DELT4  =  DELTA 
C ( 1 )  =  HS ( 1  *  1  ) 

C  (  2  )  =  C  (  1 ) 

C  (  3  )  =  H  (  2  »  2  ) 

C  (  4  )  =  H  (  3  *  2  ) 

C ( 5 )  =  C.O 
GO  TO  122 
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SUBROUTINE  SOLVE  ...(CONT'D) 

77  IF (X(J)“X{JL) -DELTA) 11 *12*13 

12  K J=K j+1 
KA J=K J 
JL  =  J 

DELT1  =  DELTA 
DELT2  =  DELTA 

11  IF(DELTA-(YY-Y(J) ) )  14  » 1 5  *  1 6 
16  IF(M- 1)20*21*22 

20  jAR1= J-l 
JAR2= JAR1 
N‘  =  vi+1 

GO  TO  40 

21  JAR1=J 
JAR2= J-l 
M  =  M+ 1 

GO  TO  40 

22  JAR1=j 
J  A  R  2  =  J 
M  =  M+l 

40  C ( 1 ) =HS ( JAR l- 1 *  2 ) 

C(2)=HS( JAR2-2>3) 

C ( 3 ) =  H ( I+1*KJ) 

C ( 4 ) *  H ( I  +  2  »  K J ) 

C ( 5 )  =  0.0 

DELT3  =  DELTA  -  ( YY  -  Y( j)  ) 

DELT4  =  DELTA 
GO  TO  122 

13  M  =  0 

IF  ( j  -  KGP ( KQA  )  ) 27  *26  »  27 

26  KAJ  =  KJ  -  1 

27  C(1)=H(I+1,kAj) 

C ( 2 ) =  H ( I  +  1*KAJ-1  ) 

C ( 3 ) =HS ( J  » 2  ) 

C ( 4 ) =  HS ( J  *  3  ) 

C  C  5 )  =  0.0 
KAj=K J 

DELT3  =  DELTA 
DELT4  =  DELTA 
IF  ( J  -  KGP ( KQA  )  )  29  >28  *29 

28  JA  =  J  -  1 
KQA  =  KQA  +  1 

29  JA  =  J 

DELT1  =  X(J)  -  X ( jA  -  1) 

DELT2  =  DELTA 
JA  =  J 

GO  TO  122 

14  YY=YY~DELTA 
1  =  1  +  1 
GO  TO  11 

122  A( 1*2)  = 
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SUBROUTINE  SOLVE  ...(CONT'D) 

A ( 1 >3 ) =  A ( 1 >2 ) #*2 
A { 2  *  2 )  =  - ( DELT 1+DELT2 ) 

A(2*3)=A(2>2)**2 
A  (  3  »  4  )  =  -DELT3 
A(3f5)aA(3#4)**2 
A ( 4  » 4 )  =  “ ( DELT3  +  DELT 4 ) 

A (4,5 )*A ( 4 >4 ) **2 

CALL  ARRAY (2>NE  » N  E  »  5  *  5 »  S  A  »  A ) 

CALL  S I MQ ( SA ♦ C  *  NE  >  0 ) 

H$ ( J » 1 ) =C (  1  ) 

100  CONTINUE 
GO  TO  18 
13  WRITE  (6*114) 

114  FORMAT ( 'O' » 1 0  X  > 'UNEXPECTED  HAPPENED  IN  SUBROUTINE 
*  SOLVE' ) 

18  RETURN 
END 
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SUBROUTINE  ARRA 


SUBROUTINE  ARRA Y ( MODE » NT > NB > NC »M > S » D ) 

DIMENSION  S(25)  *  D ( 5  *  5  > 

COMMON  H(30»60)  ♦  HS ( 60  »  3 )  » X  t  60  )  * V ( 60 )  »  RES ( 50Q )  »CHK(60 
*•160)  »  KGP  (  2  *5  )  , 

1NDP ( 60 ) 

COMMON  N  >Q » DELT A 
120  IP=0 

NM  =  0 

DO  130  K  = 1  * N B 
DO  130  L  =  1*5 
IP*IP+1 

S( IP)  =  D(L»<) 

130  CONTINUE 
140  RETURN 
END 
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SUBROUTINE  SIMQ 


SUBROUT  I N  E  S I  MG ( A  » B  # NZ *  KS ) 

DIMENSION  A (25)  *  P ( 5 ) 

COMMON  H  (  3  0  »  6  0  )  >I~S  (  60  »3  )  » X ( 60  )  > Y  (  60 ) > RES ( 5  00 )  »CHK ( 60 
*•160)  >  KGP (25)  , 

1NPP ( 60 ) 

common  n*q» delta 

TOL  =  0.0 
KS  =  0 
JC=-NZ 
D065 JD= 1 »NZ 
J Y= JD+1 
JC= JC+NZ+ 1 
RIGA  =  0. 

IT  =  JC  -  JD 
D030  I LT  =  JD » NZ 
10=  IT  +  ILT 

IF(A?S(RIGA)-ABS(A( 10 )  )  ) 20  *  30  » 30 
2  0  RIGA  = A ( I  0 ) 

IMAX  =  ILT 
30  CONTINUE 

IF(ABS(BIGA)  -  TOL ) 35  *  3  5  » 40 
35  KS  =  1 
RETURN 

40  1 1  =  JD+NZ* ( JD-2 ) 

IT  =  IMAX  -JD 
DO 50  K  =  J  D  »  N  Z 

1 1  =  1 1+NZ 

12  =  II  +  IT 

SAVE  =  A (  1 1 )  , 

A  (  1 1 )  =  A  (  I  2  ) 

A (12)  =  SAVE 
50  A(  ID  =  A(  ID  /  RIGA 
SAVE  =  B (  I  MAX ) 

B ( IMAX)  =B(JD) 

B(JD)  =  SAVE/R IGA 
I  F ( jD-NZ ) 5  5  »  70  *  5  5 
55  IQS  =  NZ*  ( JD-1 ) 

D065 I X* JY  *  NZ 
IXJ  =  IQS  +  IX 
I T  =  JD- 1 X 
D060  JX  =  JY  » NZ 
IXJX  *NZ* ( JX-1 )  + 1  X 
JJX  =  IXJX  +  IT 

60  A (IXJX)  =  A (IXJX)  -  (A(IXJ)*A( JJX) ) 

65  » (  IX) 30 ( IX )-  (B( JD)*A( IXJ) ) 

70  NY=MZ-1 
IT=NZ*NZ 
D080  JD  =  1 »  NY 
IA=IT-JD 
IB=NZ-JD 
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SUBROUTINE  SIM0  .  «.(CONT»D) 


I  C  =  N2 

DO  SO  K  =  1 »  JD 

B { 19)  =  B (  IB)  -  A ( I A ) *3 ( IC) 
IA  =  I A-NZ 
80  IC  =  IC  -  1 

RETURN 

END 


APPENDIX  B 


SAMPLE  CALCULATION  FOR  DISCHARGE 


Sample  Calculation  for  Discharge 


Applying  Darcy's  law  for  the  calculation  of  seepage 
discharge,  per  unit  length  in  an  isotropic  material  caused 
by  a  difference  in  head  H,  and  using  the  notations  shown 
in  Fig.  2.1,  it  foil  ows  : 


AQ 


(1) 


q  =  Np  AQ  and  H  =  Ah 


where  Nc  and  N.  are  the  number  of  flow  channels  and  head 
F  d 

drops  respectively. 

From  equations  (1)  and  (2) 


q 


K.H. 


An 

As 


For  a  squared  flow-net  =  1.  Therefore, 


q  =  K.Ah  x  Np  C3) 

where  K  =  coefficient  of  permeability 

Ah  =  head  drop  for  each  eq ui po tent i al  line. 


« 


B-2 


The  discharge  calculations  are  given  here  for  the 

case  correspond!' ng  to  Fig.  3.1  and  q  is  denoted  by  q 

M  J  Mmax 

Referring  to  Fig.  3.1,  it  follows 

H  =  1000 
d  =  0 
Ah  =  50 
=  0.05H 
=  2.635. 


There  fore , 


qmax  =  k  x  °-05H  x  «F 

=  k  x  0.05H  x  2.635 


or , 


q 

Mmax  _ 
KH 


2.635  x  0.05 


=  0.1317. 


This  is  the  dimensionless  representation  of  the  seepage 
discharge. 


